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1 . Introduction 

This report compiles the various research activities conducted under the auspices of the 
NASA grant NAG3-1026, "Numerical Investigation of Hot Gas Ingestion by STOVL Aircraft" 
during the period of April 1989 to April 1996. The effort involved the development of multigrid 
based algorithms and computer programs for the calculation of the flow and temperature fields 
generated by Short Take-off and Vertical Landing (STOVL) aircraft while hovering in ground 
proximity. Of particular importance has been the interaction of the exhaust jets with the head wind 
which gives rise to the hot gas ingestion process. The objective of new STOVL designs is to 
reduce the temperature of the gases ingested into the engine. 

The present work describes a solution algorithm for the multi-dimensional elliptic partial 
differential equations governing fluid flow and heat transfer in general curvilinear coordinates. The 
solution algorithm is based on the multigrid technique (Brandt [1,21, Vanka [3], Vanka et al. [4] 
and Sockol [5]) which obtains rapid convergence of the iterative numerical procedure for the 
discrete equations. Initial efforts were concerned with the solution of the Cartesian form of the 
equations. This algorithm was applied to a simulated STOVL configuration in rectangular 
coordinates (Flicker, Holdeman and Vanka [15]). In the next phase of the work, a computer code 
for general curv ilin ear coordinates was constructed. This was applied to model STOVL geometries 
on curvilinear grids. The code was also validated in model problems. In all these efforts, the 
standard k- e model (Launder and Spalding [6])was used. 

2. Problem Considered 

The ingestion of hot gases into the engine inlets is an important consideration in the design 
and operation of Short Take-off and Vertical Landing (STOVL) aircraft hovering in ground 
proximity. Such an ingestion can cause significant problems for the engine performance including 
reduced thrust and compressor stalls. These problems involve many hazards for the pilots 
including very hard landings. 

A number of interesting fluid dynamic effects have been identified when the lift jets of 
STOVL aircraft impinge on a ground surface. First, as the lift jets expand, they entrain the 
surrounding fluid, causing a negative pressure underneath the fuselage and a loss in lift. As the jets 
impinge on the ground and spread radially outwards, the exhaust jets further entrain external fluid 
and increase the loss in lift. In a multiple-jet aircraft such as the McDonnell Aircraft Company 
Model 279-3C, these exhaust jets collide with each other and turn upwards to form an up wash 
foun tain , as shown in Figure 1. This fountain flow has two major effects on the STOVL aircraft 


1 



dynamics. First, an increase in lift force is caused when the fountain impinges on the aircraft 
fuselage. The recovery in lift is a positive effect of the up wash flow. However, this impinging 
fluid can also flow along the fuselage surface and eventually make its way into the engine inlets, 
which is the primary source of hot gas ingestion. Because the temperature of the fountain flow is 
much greater than that of the ambient air, its ingestion by the engine can reduce the power and 
cause thermal stresses in engine components. In addition to the fountain flow. Figure 1 also shows 
another mechanism for the hot gas ingestion resulting from the interaction of the forward-moving 
exhaust jet with the head wind. When the head wind and the exhaust jet collide, a stagnation region 
is formed, and the exhaust jet is turned into a ground vortex. This ground vortex can subsequently 
be re- ingested by the engine, resulting in a further loss in power. Far field hot gas ingestion is 
characterized by a strong dependence on head wind velocity and aircraft height. It can be 
recognized that the overall flow field governing these fluid dynamic interactions is strongly 
three-dimensional, as well as turbulent 

The objective of this study is to develop a computational capability to investigate flow Fields 
generated by the impinging lift jets in the presence of cross flow, as applied to the STO VL flow 
geometry. Numerical simulation of the complete STOVL aircraft flow field requires the solution of 
nonlinear partial differential equations that govern the transport of mass, momentum, and heat in 
three dimensional space with boundary conditions that describe the lift jets, the aircraft geometry, 
and the head wind. 

In reality, the flows of practical relevance are almost always turbulent; this means that the 
fluid motion is highly unsteady and three-dimensional. Due to these complexities, the turbulent 
motion and the heat and mass-transfer phenomena associated with it are extremely difficult to 
describe and thus to predict theoretically. Yet, one of the basic tasks of fluid engineering is to 
predict fluid flow phenomena. Due to the fact that "predictions" by way of experiments are usually 
very expensive, calculation methods are in demand. 

A complete simulation of a turbulent flow requires the solution of the time dependent 
Navier- Stokes equations. Such a simulation is referred to as a direct numerical simulation (DNS) 
(Orszag and Patterson [7]). A direct numerical simulation is capable of resolving the temporal and 
spatial dependencies of the turbulent flow field. However, the ability to resolve the significant time 
and length scales is achieved at the cost of computational effort As an example, a fully turbulent 
isothermal channel flow calculation using a DNS required 250 hours of CPU time on a Cray X-MP 
super computer (Kim et al. [8]). At present time, DNS is not a viable option for the solution of 
practical flow fields. 
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The solution of time averaged equations is currently the most widely used means of 
simulating an engineering flow problem [6], In solving a set of averaged equations, a temporally 
and spatially resolved solution is replaced by a stochastic description of the flow. The dependent 
variables of the governing equations are time averaged quantities. 

The solution of averaged equations is obtained only after certain closure hypotheses have 
been invoked. The closure hypotheses of turbulence models relate the correlation between the mean 
fluctuations (Reynolds stresses) to the time averaged flow variables. The degree of sophistication 
of this closure is dependent on the flow fields which are to be studied. A review of the various 
models which are in use for predicting internal flows is given by Nallasamy [9], For wall bounded 
flows, Patel et al. [10] give an analysis of various low Reynolds number k- e turbulence models. 
A review of Reynolds stress models is given by Speziale [1 1]. The standard k-e turbulence model 
is used in the current study. 

3 . Objectives of Current Research 

The scope of the current research project was to develop efficient computational techniques 
for elliptic fluid flow equations and apply them for studying the mechanisms and extent of hot gas 
ingestion by Short-Take-off and Vertical Landing (STOVL) aircraft. The efficient computational 
technique was based on the solution of the Reynolds- averaged equations by a multignd technique. 
In a multigrid procedure, in contrast with a single grid procedure, several levels of fine and coarse 
grids are used to discretize the partial-differential equations. The solution of the discrete equations 
is sought on the finest grid, but use is made of the coarse grids to accelerate the solution process. 
The coarse grids are also used to progressively generate better initial guess for the finer grid. 

The general objective of the present research was to develop a computer code capable of 
simulating flow in complex geometries with embedded blocked regions (obstacles and baffles) as 
well as arbitrary locations where mass and momentum can be injected. These features are 
necessary to simulate the STOVL configuration and the dynamic interactions occurring during the 
hovering stage. The construction of the eventual code included an intermediate stage of considering 
a rectangular (Cartesian) geometry and validation of the algorithm in model problems. 
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4 . Present Contribution 

During the course of this research, several versions of the computational algorithm were 
developed and tested. Initially, a version limited to rectangular coordinates was developed and 
applied to compute twin-jet impingement (Pegues and Vanka [12]) and flow over a model STOVL 
aircraft (Tafti and Vanka [13,14]). This version was subsequently used by Fricker, Holdeman and 
Vanka [15] to study the effects of geometrical and flow parameters on the hot gas ingestion. 

The curvilinear code development consisted of first developing an algorithm based on 
collocated arrangement of the variables. The collocated arrangement has several advantages over 
other formulations of curvilinear grid equations. The multigrid technique was developed for 
collocated curvilinear grids as a part of this study. The code and the algorithm were validated in a 
number of model problems as well as in a STOVL configuration (Smith [16], Smith and Vanka 
[17], and Smith et al. [18]). However, during further trials with this code in more realistic STOVL 
configurations, difficulties related to convergence at the end of the calculation were encountered. 
This was attributed to the complex nature of the multigrid cycling procedure on collocated grids. 

As a result, a new staggered arrangement based on previous work by Maliska and Raithby 
[19] was considered. This scheme was subsequently developed for a multigrid sequence and 
validated in complex flows (Cope et al. [20], Wang et al. [21] and Wang [22]). Based on the 
success of these calculations, the STOVL calculations were repeated with this method. However, 
similar difficulties as encountered with the collocated arrangement were also encountered with this 
arrangement. 

The STOVL problem has several complexities from computational viewpoint that were 
probably the reason why good performance was not observed in this flow while expected fast rates 
of convergence were observed in other fairly complex flows. The specific issues of STOVL are: 

a) large cell aspect ratios are necessary to represent the large stream wise extent of the flow domain; 

b) large rates of strain and turbulence production are encountered at the impingement sites of the 
jets resulting in questionable performance of the wall-functions; 

c) large pressure build-up at stagnation regions (of jets) resulting in severe pressure gradients at the 
bottom boundary and inconsistencies from pressure extrapolation to the walls; 
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d) large-scale unsteady behavior of the jets in the impingement region causing difficulties due to 
representation as a steady flow. 

Further research is needed to address these issues and to identify the precise cause for these 
convergence difficulties. To this end, methods based on unstructured triangular grids may be more 
beneficial. Also, unsteady calculations in the context of Very Large Eddy Simulations (VLES) may 
appear to be promising for future studies. 
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& Exhibit, January 6-9, Reno, Nevada, AIAA Journal of Aircraft, vol. 31,1, pp. 236-242. 
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Geometries”, AIAA-92-0096, 30th Aerospace Sciences Meeting & Exhibit, January 6-9, 
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NUMERICAL STUDY OF TWIN-JET IMPINGEMENT 
UPWASH FLOW 


W. J. Pegues and S. P. Vanka 

Department of Mechanical and Industrial Engineering 
University of Illinois at Urbana-Champaign 
Urbana, Illinois 


Abstract 

Two horizontally spaced jets impinging normally on a flat surface 
create a fountain upwash flow due to the collision of the radially 
flowing wall jets. This fountain flow is of importance to the dynamics 
and propulsion of Short Take-off and Vertical Landing(STOVL) 
aircraft . The fountain flow influences the lift forces on the aircraft and 
the ingestion of hot gases and debris by the engine inlet In this paper, 
a multigrid based finite-difference numerical procedure has been 
applied to solve the equations governing this three-dimensional flow. 
Tlie calculations have been performed using a reasonably fine finite- 
difference mesh and the results have been compared with experimental 

data of Saripalli (1985). The standard k-e turbulence model has been 
used. Comparisons with experimental data reveal that while the mean 
velocities are predicted with reasonable accuracy, the turbulent kinetic 
energies are seriously in error. The reasons for this discrepancy could 
be the intense unsteadiness and large-scale structures of the flow in the 
near wall region, which can not be captured well by any Reynolds- 
averaged turbulence model. 

Introduction 

For the past three decades, there has been a significant amount of 
research into the understanding and predictability of the flow fields 
that govern the dynamics of Vertical Take-off and Landing (VTOL) 
aircraft (Mitchell, 1985 and Gilbert, 1984). The impingement of the 
lift jets on the pound produce a complex three-dimensional flow 
consisting of radial wall jets and an upwash flow. The upwash flow 
subsequently impinges on the fuselage and impacts the lift and 
moments on the aircraft Also, part of tire hot gases impinging on the 
fuselage make their way into the engine inlet and can cause thermal 
loading on the engine components (Kuhn, 1982; Mitchell, 1985 and 
VanOverbeke and Holdeman, 1988). Therefore, detailed studies of the 
fountain upwash flow are essential to avoid adverse effects due to the 
impinging lift jets. Also, the upwash flow presents itself as an 
interesting fluid flow problem from the viewpoint of basic fluid 
mechanics. 

In recent years, several experimental and computational studies have 
been conducted to study the fluid mechanics of single and twin 
impinging jets. Agarwal (1983), Kuhn and Eshelman (1985) and 
Kotansky (1983) have provided reviews of current understanding and 
computation of these flow fields. Several complexities of the flow 
have been pointed out. Foremost of these is the three-dimensionality 
of the flow which makes numerical computations difficult and time 


consuming even on today's supercomputers. Also, flow visualization 
experiments(Kibens eL al.,1987) indicate that the jets and the fountain 
region are highly unsteady and are dominated by large scale turbulent 
structures. Thus, attempts to model the effects of turbulence through a 
Reynolds-averaged approach can be inaccurate and frustrating. In 
earlier works of Agarwal and Bower(1982) and Kotansky and 
Bower(1978), the flow was simplified to be two-dimensional and the 
governing equations were solved by a finite-difference calculation 
procedure. Turbulence models based on solution of partial-differential 
equations for turbulence kincric -energy and another variable such as 
the turbulence dissipation rate have been used to represent the effects 
of turbulence. Similar computational studies have also been conducted 
by other researchers (e.g. Hwang et aL, 1990, and Mitchell, 1985). 
Despite the difficulties in accurately modelling turbulence, the 
Reynolds-averaged approach continues to be the preferred approach 
because time-accurate calculations of the unsteady flow are 
significantly more expensive (Childs and Nixon, 1986; Rizk and 
Menon, 1988 and Childs, 1990). 

A primary difficulty in performing three-dimensional fluid flow 
calculations (Bower, et al, 1979) is the inability to resolve the flow 
domain with enough mesh points such that the numerical errors in the 
solution are essentially negligible. In the present work, we have used 
a multigrid based calculation procedure to compute the three- 
dimensional flow field generated by the impingement of two 
horizontally-spaced jets. The multigrid technique possesses the 
advantage that the discretized elliptic partial differential equations can 
be solved in a relatively small number of iterations even when the 
mesh is significantly refined. Therefore, calculations with fine mesh 
resolutions can be performed in relatively small computer times 
leading to the possibility of an error ffee solution.The multigrid 
concept is explained in several earlier works (Brandt, 1984 and 
Hackbusch, 1985), and our implementation is documented in Vanka 
(1988), Vanka eud. (1988) and Claus and Vanka (1990). 

Flow Configuration 

Recently, Saripalli and Kroutil (1985) and Saripalli (1981,1985) 
conducted detailed experiments on the twin jet impingement flow. 
Initially, flow visualization experiments were conducted in a closed 
rectangular tank with two water jets attached to a baffle plate (Fig! 1). 
The water jets impinged normally on the bottom of the tank and 
created a fountain upwash flow. Subsequently, a laser velocimeter 
was used to map the mean and fluctuating velocity components. The 
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effect of jet spacing and the height of the jets from the ground plane 
were also studied. Although the experiments were isothermal, these 
data provide valuable information towards developing and validating 
turbulence models for engineering calculation of STOVL aircraft flow 
processes under non isothermal conditions. 



Equations and Solution Procedure 

The equations governing a three-dimensional steady, elliptic flow of 
constant density and viscosity can be written as 

dui/3xi = 0 (1) 

3/9xj (ui uj) = - dp /p9xi + v 9y9xj (chii/ftxj) , i =1,3 (2) 

where summation is implied by the repeated indices, ui is the 

component of velocity in the xi direction and p is the pressure, v is the 
molecular viscosity. In the Reynolds-averaged approach, the velocities 
are decomposed into mean and fluctuating components and a 
turbulence model is used to repres ent the correlations of the fluctuating 

velocities. The most popular of these models is the k-e model and a 
large amount of research effort has been devoted to its development 
and validation. Higher order models based on stress transport 
equations (Launder, 1989) are also available, but at this stage further 
development efforts are needed before they are applied to three- 
dimensional flows such as the one considered here. 

In the k-e model, two other partial-differential equations are solved 
representing the transport of turbulence kinetic energy and its 
dissipation rate. The model assumes that the turbulent stresses can be 
related to the mean strain rates through a turbulent viscosity. This 
turbulent viscosity is evaluated from the distributions of turbulence 
kinetic energy and its dissipation rate. Validation of this model in a 
wide variety of flows has been summarized by Rodi (1984), amongst 
others. Although there are serious doubts on the validity of the 
assumptions made in deriving this model, it currently represents the 
best compromise for engineering flows. 


Computational Details 

The equations are discretized by a first/second order accurate hybrid 
differencing scheme (Spalding, 1972). Although higher order 
schemes can be used to improve the accuracy of the discretization, 
such schemes can be non conservative and may not provide the 
expected accuracy for complex recirculating flows (Vanka, 1987). In 
the current differencing scheme, the convective fluxes are discretized 
such that at high Peciet numbers (>2), the fust order scheme is used, 
whereas for Peciet numbers < 2, the second order accurate central 
differencing scheme is used. This hybrid scheme currently appears to 
be the only robust scheme for incompressible, recirculating flows as 
other proposed schemes have been observed to cause difficulties in 
convergence and overshoots in the flow variables (Syed et. al., 1985). 

The discrete equations arc solved by a multigrid solution algorithm. A 
fine grid consisting of 64 x 32 x 32 finite -difference cells is used to 
represent the important flow domain. In the experiments of Saripalli 
(1985), the rectangular tank is closed at all the ends, except for a small 
outlet to remove the jet fluid. We represented this geometry by placing 
a baffle inside a rectangular solution domain and included the jets as 
mass injections to internal cells below the baffle. The jet fluid was 
extracted above the baffle as a mass and momentum sink. Only one- 
fourth of the total geometry was simulated because of symmetry about 
the span wise mid-plane and the mid-plane of the jets.The jets were 
considered to be rectangular because it was not currently possible to 
represent mass injection from a circular orifice. 

The multigrid calculation consisted of four grid levels starting from the 
coarsest grid of 8 x 4 x 4 cells. This grid was progressively doubled 
(mesh spacing is halved in each direction) to a 64 x 32 x 32 grid. The 
two jets were spaced nine jet diameters apart, corresponding to the 
first of the three tests of Saripalli. The jet was discretized with sixteen 
cells in the x-direcrion and eight cells in the z-direction. The fountain 
region between the jets was represented by thirty two cells in the x- 
direction. In the y-direction, a non-uniform grid with closer spacing in 
the wall region was used. The calculations were iterated until the sum 
of the absolute local mass residuals was less than one percent of the 
jet mass flow. 

The multigrid procedure used here is identical to the one reported 
earlier (Vanka et. al, 1988). The inner relaxation of the continuity and 
momentum equations is performed simultaneously by a coupled 
Gauss-Seidel iterative procedure. The turbulence equations are solved 
decoupled from the momentum and continuity equations and interact 
with the flow equations through the turbulent viscosity. A Full 
Multigrid (FMG) algorithm is used in which the solution is initiated at 
the coarsest grid and progressively developed for finer grids. The 
nonlinear version of the multigrid procedure is used on any given grid 
level, with interpolations between grids dictated by a V-cycle. Details 
of these features are given in Vanka et al., 1989. 

Results 

Figure 2 shows the observed convergence history of the mass 
residuals for S/D = 9.0 and H/D = 3.0 (H = height of jets; S = spacing 
of jets; and D = jet diameter). It can be seen that the procedure 
converges rapidly to a satisfactory solution, requiring approximately 
80 iteradoos to reach the one percent level from the initial level of 300 
percent. As a result of this superior rate of convergence, the 
calculations were completed in approximately 8 minutes on a CRAY-2 
computer. Although the grid employed here has not been proven to 
provide error- free results, it is sufficient for conducting parametric 
studies as well as for any future studies with higher-order 
discretization schemes. Also, calculations with another finer grid level 
will be attempted in the near future. 

A systematic comparison of the present calculations with data of 
Saripalli (1985) has been performed both in the jet and in the fountain 
regions. Saripalli measured the time-averaged mean velocities in the 
axial and radial directions of the jet and the fountain (his notations 
were interchanged in the jet and fountain regions). Also measured 


13 





Fig. 2 Convergence History for Sum of Mass Residuals 


were the fluctuating components ( u’ 2 and v 2 ) from which a 
turbulence kinetic energy could be derived (assuming that w'2 is 
equal to the average of u' 2 and v r 2 ). The turbulence kinetic energy 
from the solution of the partial differential equation was compared to 
this value. 

Figure 3 provides a comparison of the profiles of axial velocity in the 
jet region at three different positions from the ground. The heights 
selected are closer to the ground because the flow docs not develop 
much in the initial region and the agreement with the experiments is 
very good in that region. For the three different heights shown in 
Figure 3, we observe that the agreement at y/D = 1.0 is fairly good. 
The minor differences in the central core region and at the edges of the 
jet can probably be reduced by more accurate prescription of the inlet 
conditions and further mesh refinement. At y/D = 0.3. the calculated 
jet velocities are lower (considering the magnitude) in the core region 
and higher at the edges. This region is affected by the impingement 
region and the errors in predicting the latter are propagated upstream. 
Closer to the wall, there are appreciable differences, particularly to the 
left of the jet (on the other side of the fountain). It is necessary to point 
out that this region is affected by the way the boundary conditions are 
prescribed for extracting the jet mass, in the current calculations, the 
flow is extracted as a mass sink in the region above the baffle. This 
makes the flow turn upwards towards the mass sink, influencing the 
region left of the jet However, in the region adjacent to the fountain, 
the calculated velocities show reasonable agreement. 

Figure 4 compares the calculated radial velocities in the wall jet region 
with the measured values. At y/D = 1.0, the radial velocities are 
relatively small, and the qualitative nature of the velocity profile is 
accurately captured by the calculations. At y/D of 0.30, the flow on 
the fountain side is overpredicted but the flow on the free stream side 
is under-predicted. However, qualitatively, the profile is well 
predicted. At y/D of 0.05, the velocities are again under-predicted, 
implying that the wall jet profile is flatter than that in the experiments. 
This may either be due to insufficient grid refinement in the wall jet 
region or due to the inaccuracies in the turbulence model. 

In Figure 5, the measured and calculated kinetic energies are 
compared. The trends observed for the velocities are also seen here, 
but the discrepancies in the near-wall region are much larger. The jet 


shear layer is fairly well captured at y/D of 1.0. However, at y/D of 
0.30, the predicted profile of kinetic energy is flatter than the 
experiments and significant differences exist. The disagreement is 
even larger for y/D of 0.05. The experiments indicate large values of 
kinetic energies on both the fountain side and the free stream side of 
the jet. However, the calculations substantially under-predict these 
values, indicating a flat distribution of turbulence energy. Near the 
wall, the jet impingement region is very violent and unsteady (Childs 
and Nixon, 1986). This produces large turbulent fluctuations, which 
are reflected in the measurements as high turbulence kinetic energies. 
Such fluctuations are not well represented by a Reynolds-averaged 
turbulence model. Much of the discrepancy observed in the 
impingement region is believed to be due to the inherently large scale 
nature of the turbulent flow. However, some improvement in the 
results may be achieved if a finer mesh is used closer to the wall. 

The velocity field in the fountain region is compared with experimental 
data, in Figures 6 and 7. Figure 6 shows that the horizontal velocity 

(u) is well predicted in the present calculations. These velocities are 
approximately 20 percent of the jet velocity, and the agreement is 
fairly good except at y/D of 0.05. For this location, the axial velocity 
close to the jet is under-predicted. The agreement in vertical velocity 

(v) , shown in Figure 7 is better than that for the horizontal velocity. 
The location x/D of zero corresponds to the center of the fountain 
where the vertical velocity is largest. This value is well predicted. 
However, significant differences exist between the calculated and 
measured turbulence kinetic energies shown in Figure 8. The 
measurements show large values of turbulence kinetic energy in the 
center of the fountain region. This turbulence is produced by the 

collision of the two wall jets. The calculations performed with the k-e 
model significantly under-estimate these energies. Significant 
differences exist at ail three values of y/D. The disagreements 
observed in turbulence energies can not be improved easily by the 

standard k-e model. Recent studies (Childs, 1990) indicate that several 
corrections which include effects of streamline curvature, anisotropy 
and large scale mixing are necessary before good predictions are 
possible. Alternatively, tune accurate simulations which solve for the 
large scale structures may provide more hope in representing the 
physics of the flow. 

In the present study, calculations for the other two sets of data of 
Saripalli (1985) in which the S/D and H/D are varied were also made 
and compared with the measured values. The agreement between the 
calculations and experiments is similar to that presented above. 
Therefore, for the sake of brevity, the other results are not presented 
here. 

Summary 

This paper demonstrates the applicability and potential of a multigrid 
procedure for efficiently solving the fountain up wash flow generated 
by the lift jets of Short Take-off and Vertical Landing(STOVL) 
aircraft. Numerical computations have been performed for the 
experimental configuration of Saripalli (1985) with a finite-difference 
grid consisting of 64 x 32 x 32 cells. Comparisons of the calculated 
velocities and turbulence energies reveal that while the velocities agree 
reasonably well with measured values, the turbulence kinetic energies 
are in significant error. Although some improvement can be achieved 
by further mesh refinement, it is believed that the fundamental cause 
for the disagreement may be the unsteadiness of the impingement 
region which can not be well represented by the Reynolds-averaged 
concept. Therefore, time-accurate simulations, based on the concept of 
simulating the large scale structures (LES) may be more appropriate 
for accurate calculations of this flow. 
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Fig. 3 Axial Velocity Across the Impinging Jet 


Fig. 4 Radial Velocity Across the Impinging Jet 
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Fig. 5 Turbulent Kinetic Energy Across the 
Impinging Jet 
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Fig. 6 Axial Velocity in the Fountain 
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Fig. 7 Radial Velocity in the Fountain 


Fig. 8 Turbulent Kinetic Energy in the Fountain 
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Abstract 

Ingestion of hot exhaust gases by the engines of Short Take-off 
and Vertical Landing (STOVL) aircraft has been an important 
research problem for several years. The hot gas environment 
around STOVL aircraft is three-dimensional and turbulent. In 
this study, the Navier-Stokes equations governing the hot gas 
ingestion flow field are solved by an efficient finite-difference 
calculation procedure. The complete geometry including the 
head wind and the fuselage is simulated. Four demonstration 
calculations with variauons in the height of the fuselage and the 
head wind velocity are presented. It is shown that the 
calculation procedure efficiently provides a soluuon to the 
governing equauons and produces realisuc descriptions of the 
flow and temperature fields. 

1 . Introduction 

A number of interesting fluid dynamic effects have been 
identified when the lift jets of Short Take-off and Vertical 
Landing(STOVL) aircraft impinge on the ground surface 
(Kuhn [1], Kotansky [2], Kuhn and Eshelman [3], Agarwal 
[4]). First, as the lift jets expand, they entrain the surrounding 
fluid causing a negative pressure underneath the fuselage and a 
loss in lift. As the jets impinge on the ground and spread 
radially outwards, the wall jets further entrain external fluid and 
increase the loss in lift. In a multiple-jet configuration, these 
wall jets collide with each other and turn upwards to form an 
upwash fountain. This fountain flow has two main effects on 
the STOVL aircraft dynamics. First, an increase in lift force is 
caused when the fountain impinges on the aircraft fuselage. The 
recovery in lift is a positive effect of the upwash flow. 
However, this impinging fluid can also flow along the fuselage 
surface and eventually make its way into the engine inlets. 
Because the temperature of the fountain flow is much hotter 
than the ambient air, its ingestion by the engine can reduce the 
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power and cause thermal stresses in the components. In 
addition to the fountain flow, another mechanism for hot gas 
ingestion results from the interaction of the forward moving 
wall jet with the headwind.When the headwind and the wall jet 
collide, a stagnation region is formed and the wall jet is turned 
into a ground vortex [3]. This ground vortex can subsequently 
be re-ingested by the engine, resulting in a further loss in 
power. However, the hot gas ingestion process depends on 
where and how the inlets to the engine are located with respect 
to the approach flow. For strong headwinds, the ground vortex 
can be pushed behind the engine inlets such that no hot gases 
are ingested. It can be recognized that the overall flow Field 
governing these fluid dynamic interactions is strongly three- 
dimensional as well as turbulent. 

In the last three decades, a number of studies addressing the 
above issues have been conducted(see [4,5] for good reviews). 
A majority of these studies were experimental in nature 
although some analytical studies of the individual processes, 
such as single jet impingement have also been conducted 
[6,7,8]. These studies have identified the fundamental 
processes and to a certain extent, quantified the effect of the 
parameters on hot gas ingestion and fountain flow. However, 
detailed measurements of the velocity and temperature fields in 
the fountain region and in the ground vortex in realistic 
operating conditions are not currently available. Recently, 
Saripalli [9,10] conducted some experiments in a water tunnel 
in which Laser-Doppler Velocimeter and Laser Imaging 
techniques were used to study the fountain region of the lift 
jets. These studies involved isothermal flows and were limited 
to the fountain formation. 

With recent advances in the development of powerful digital 
computers, it has now become possible to numerically solve the 
equations that govern the transport of mass, momentum and 
enthalpy. However, because the flow is turbulent, such studies 
must necessarily make assumptions on the macroscopic 
behavior of turbulence. Despite the inherent limitations of such 
an approach in representing all the effects of turbulence with 
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precision, it has become a useful engineering tool. Numerical 
computation of the complete STOVL aircraft flow field requires 
the solution of nonlinear partial differential equations that 
govern the transport of mass, momentum and heat in three- 
dimensional space with boundary conditions that describe the 
lift jets, aircraft fuselage and the head wind. Because of the 
complexity of such solutions and limited computer capabilities, 
until recently most numerical studies were limited to two 
dimensions in which simplifications in the spanwise direction 
were made. Kotansky and Bower [11] were probably the first 
investigators to apply a Navier-Stokes analysis to the study of 
planar jet impingement of relevance to VTOL aircraft. In their 
analysis, a one equation model for turbulence kinetic energy 
was integrated with equations describing the stream function 
and transport of voracity. A numerical solution of an impinging 
jet was obtained and compared with experimental data. This 
study was followed by Agarwal and Bower [12], who 
improved the turbulence model by considering an additional 
equation for the turbulence dissipation. Because the improved 
model explicitly calculated the local length scale of turbulence, 
the flowfields predicted by their model displayed better 
agreement with experimental data. Recently, other studies have 
appeared which also have numerically solved the equations for 
single and twin jet impingement flows using finite-difference 
techniques [e.g. 13,14,15] 

Although the hot gas ingestion process results primarily from 
the impingment of the lift jets, the eventual ingestion is dictated 
by the complete flow dynamics including the formation of the 
fountain up wash and its interaction with the headwind. For 
stronger headwind speeds, it is possible to push the ground 
vortex downstream of the engine inlet and thus completely 
avoid its contribution to hot gas ingestion. Also, the collision of 
the wall jets and the fountain formation can be affected by the 
headwind and the geometry of the fuselage under surface. In 
addition, placement of Lift Increasing Devices (LIDs) [3] and 
other obstructions can divert the fountain flow from the engine 
inlets, thereby reducing the severity of the hot gas ingestion. In 
order to better understand all the flow processes, a complete 
Navier Stokes analysis including the far field condition of a 
prescribed headwind is necessary. Although such an analysis is 
computationally very intensive, the benefits are significant for 
evaluating the effects of geometric and flow parameters. 
Recently, an important contribution has been made in this 
direction by VanOverbeke and Holdeman [16]. In this study, 
the complete fuselage, headwind and multiple impinging jets 
were numerically simulated and the temperature fields close to 
the engine inlet face were studied. By varying the head wind 
and the impingement height, they were able to alter the 
ingestion pattern and quantify their effects. A finite -difference 


computer program [17] was used and the flow domain was 
discretized into a relatively large number of finite volumes. 

A primary drawback of numerical computations of three 
dimensional elliptic fluid flows is that they can require 
substantial amounts of computer time to obtain a converged 
solution. These large computational costs often discourage 
systematic and thorough investigations of the effects of various 
parameters. However, with research into faster converging 
numerical algorithms, it is possible to substantially reduce these 
computer times. In this paper, we demonstrate the applicability 
of one such algorithm [18] which is based on the concept of 
using multiple levels of grids [19] to obtain faster convergence. 
In this study, the geometrical configuration studied by Van 
Overbeke and Holdeman [16] is considered and a set of four 
calculations have been performed. To compare our results with 
those of VanOverbeke and Holdeman [16], we have simulated 
exactly the same configurations, albeit with a finer grid, and 
compared the computer times and the flow fields. The four 
calculations in this study differ in the height of the jets from the 
ground and the ratio of the head wind to the jet velocity. 

The following sections briefly describe the solution procedure 
and present the results of the calculations. In Part 1 of this 
study, a separate experimental investigation is presented by 
McLean et. ai [20] in which the concentration of the exhaust 
jets is quantified through a Laser Imaging technique. 

2 . Governing Equations and Solution Algorithm 

The flow and temperature fields around the STOVL aircarft are 
governed by the following set of elliptic partial-differential 
equations. In the present analysis, the flow is considered to be 
steady in the time-averaged sense, and the effects of turbulence 
are captured by a turbulence model [21]. The governing 
equations can be written as follows: 

Mass continuity 

^ ( p u) + ^ ( p v) + b (pw) = m (1) 

x-momentum 

d(puu) 3(puv) 3(puw) 3P 3 /_ 3 u\ 
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y- momentum 
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z-momentum 
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Turbulence energy 
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Dissipation of turbulence energy 

Jr ( P UE)+ Jr(P ve > + ^(P we ) = 


Jr( r ^) + Jr( r E 5p) 


3 3e\ „ _ e _ e 2 - 

+ 3^[ r £ 5iJ + ClG k' C 2 P¥ + mein j (6) 


Enthalpy 
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Equation of state 
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P = 


RT 


(8) 


u, V, w, are the velocities in the x, y and z 
directions,respectively. P is the pressure and R is the gas 
constant and G is the production of turbulent kinetic energy. In 
all the above equations, the diffusive coefficients contain the 
effects of laminar as well as turbulent diffusion. Standard 
values of the constants in the turbulence model [21] are used. 

The above equations are solved on a rectangular domain that 
envelopes the aircraft fuselage and the exhaust jets. Because the 
current computer program is limited to the use of rectangular 
grids, the fuselage of the aircraft was modeled to be of 
rectangular shape. However, efforts are currently under way to 
extend the program to handle grids of arbitrary inclinations 
(boundary-fitted grids [22]) which will then permit a more 
realistic representation of the aircraft shape and angles of attack. 
The boundary conditions and solution domain considered in 
current calculations arc described in section 3.1. 

The above equations arc finite-differenced by a hybrid scheme 
[23] that combines first and second order differencing for the 
total convective and diffusive operator. At high cell Peclet 
numbers, the finite-differencing scheme becomes first order 
accurate to preserve stability of the solution procedure. The 
multigrid solution technique [19] differs from the singe grid 
technique in the manner the finite-difference equations are 
solved to the required accuracy. It is well known that for elliptic 
equations, single grid techniques converge poorly when the 
finite-difference grid contains a large number of mesh points. 
This results from the low frequency errors that are slow to 
converge on a given fine grid. In the multigrid concept, a scries 
of fine and coarse grids is used and the solution is switched 
between the coarse and fine grids such that errors of all 
frequencies converge at the same rate. This novel concept was 
originally proposed for model elliptic equations, but recently it 
has become popular in solving practical fluid flows [24]. 

The principle behind the multigrid procedure may be explained 
as follows. Consider an elliptic equation that is discretized by 
finite-difference or finite element methods. Generally, an 
iterative procedure that is designed to obtain a solution to the 
discrete equations will converge rapidly in the first few 
iterations but subsequently the convergence deteriorates. It can 
be shown through formal analyses that the cause for this poor 
convergence is the presence of low frequency errors that are not 
resolved well on any given grid. However, these errors can 
converge faster if the grid is coarse. Therefore, in the muitigrid 
procedure, these errors arc interpolated to a coarser grid and 
solved. The results of this coarse grid solution arc then put 
back on the finer grid solution through extrapolation. In the 
case of a calculation with many mesh points, several layers of 
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coarse meshes can be used to obtain rapid convergence. 
However, the coarsest grid is often dictated by the constraints 
of geometry and boundary conditions. 

In the present solution scheme, we have combined the multigrid 
technique with a coupled solution scheme for the momentum 
and continuity equations. The velocities and pressures are 
obtained by a block-implicit solution of the momentum and 
continuity equations. As this procedure is explained in previous 
works by one of the authors [18,25,26], it is not elaborated 
here. However, in order to simulate the fuselage and the lift 
jets, the procedure has been extended to handle internal 
obstacles and mass and momentum injections. These features 
have been incorporated in a manner that maintains consistency 
between the coarse and fine grids, a necessary feature for the 
algorithm to eventually converge to the correct solution. It must 
be pointed out that although the coarse grids are used, the final 
solution is obtained on the finest grid in the system. Thus the 
intermediate grids are used only for accelerating the 
convergence rate and they do not influence the final converged 
solution. 

The next section describes the computational domain and a 
typical convergence history of the multignd solution procedure. 
Simulated flow and temperature fields are presented along with 
a discussion on the mechanism of hot gas ingestion for a four 
jet configuration. These results compare qualitatively with the 
flow imaging data of McLean, Sullivan and Murthy [20], as 
well as with calculations presented by Van Overbeke and 
Holdeman [16]. 

3 . Results 

3.1 Computational details and Convergence 
Histories 

In this study, four test cases have been considered, which are 
similar to those studied in ref. 16. The important flow and 
geometrical parameters are outlined in Table L The four cases 
differ in the height of the fuselage from the ground plane (h/d 
where d is the diameter of the lift jets) and the headwind to jet 
velocity ratio (Uoo/Uj). The solution domain is a three- 
dimensional wind tunnel with overall dimensions of 136d in the 
streamwise direction (x), 27 or 29d in the cross-stream 
direction (y) and 40d in the spanwise direction (z). These large 
dimensions are necessary to accurately resolve the flowfleld 
around the jets and the fuselage without any interference from 
the simulated boundaries. Figure 1 shows the computational 
grid used and the flow geometry. The geometry simulates four 


Table 1. Flow and Geometry Parameters. 


CASE 

U-/Uj 

h/d 

Maj 

Re- 

xIO- 4 

Rej 

xIO' 5 

1 

0.03 

4.0 

0.47 

3.3 

3.3 

2 

0.03 

2.0 

0.47 

3.3 

3.3 

3 

0.09 

4.0 

0.47 

9.9 

3.3 

4 

0.09 

2.0 

0.47 

9.9 

| 

3.3 


jets (due to symmetry, only half the domain is considered) 
issuing from the engine with the engine inlet located 10 jet 
diameters upstream of the fore jet axis. The distance between 
the fore and aft jet axis is 6 diameters while the spanwise 
distance between the jets is 3.2 diameters. The fuselage extends 
throughout the calculation domain and is 2.5 diameters in 
height and 2.2 diameters in width. The inlet, which is 67 jet 
diameters from the entrance of the test section, has a height of 
2.5 diameters and a width of 1.4 diameters. It is adjacent to the 
fuselage in the spanwise direction and extends 9.5 jet diameters 
up to the fore jet after which it becomes part of the fuselage. 
The calc ulati on domain is divided into 100 computational cells 
in the x direction, 44 in the y, and 48 in the z direction, giving a 
total of about 211,000 cells. The particular grid distribution 
was chosen to provide good resolution in the near field of the 
jets and at the same time keep the grid expansion ratio as small 
as possible. The grid distribution was also influenced by the 
necessity to accomodate internal obstacles and baffles such that 
grid lines bounded the surfaces of these structures. 

At the inlet of the test section, uniform velocity and temperature 
conditions are prescribed for the incoming headwind and the 
turbulent kinetic energy was assumed to be one percent of the 
mean kinetic energy. At the exit, outflow boundary conditions 
are prescribed. In the y direction, wall boundary conditions are 
prescribed. The bottom wall is assumed to have zero heat loss 
while the lop wail is assumed to be at a uniform temperature. In 
the spanwise direction, symmetry boundary conditions are 
applied on one side while an adiabatic wall is prescribed on the 
other side. Within the calculation domain, additional conditions 
are prescribed for the engine inlet and the jets. The inlet acts as 
a mass sink while the jets as mass sources. In prescribing these 
conditions it is important to conserve mass. The jet velocity 
profile is assumed to be uniform with the turbulent kinetic 
energy equal to one percent of the mean kinetic energy. 
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In the calculation procedure, the prescribed fine grid was 
coarsened to two additional levels* The calculation was started 
on the coarsest grid from prescribed initial conditions but on the 
finer grids, the starting solutions were obtained by successively 
extrapolating converged flow fields from the coarser grids. 
This Full Multigrid Cycle (FMG) provided realistic starting 
values for the calculations on the fine grids. Figure 2 shows a 
typical convergence history for.the calculation of Case 1. The 
convergence is monitored by the normalized sum of the 
absolute mass and momentum residuals over the whole solution 
domain. It is seen that the solution converges to the required 
tolerance of one percent error in approximately 150 iterations. 
In the previous study of VanOverbeke and Holdeman [16], 
which used a single grid procedure, the number of iterations 
were approximately 2000 for a 134,000 grid. This represents 
a significant reduction in computational effort. The present 
calculation required approximately 35 minutes of CPU time on 
a CRAY-2; typically this translates to a few hours of interactive 
real clock time. 

3.2 Calculated Velocity and Temperature Fields 

Hot gas ingestion is a serious problem for STOVL aircraft in 
ground proximity. Two important factors which govern hot gas 
ingestion are the proximity of the lift jets to the ground (h/d 
ratio) and the ratio of the freestream velocity to the jet velocity 
(Uoo/Uj). Another important factor, which is not considered in 
the present study, is the location of the inlet with respect to the 
lift jets and the use of shields and baffles to deflect flow of hot 
gases away from the engine inlet. When the lift jets impinge on 
the ground they stagnate and form radially flowing wall jets. 
When two opposing wall jets meet they manifest themselves 
into a fountain flow. The flow Field caused by a rectangular 
four jet configuration is governed by a strong fountain flow 
between the fore and aft jets (streamwise fountain) and also 
between the symmetrical spanwise jets (spanwise fountain). 
The streamwise fountain upwash when obstructed by the 
fuselage is forced to flow laterally outwards and up while the 
spanwise fountain upwash tends to flow under and along the 
length of the fuselage in the streamwise direction. Both these 
mechanisms contribute to the direct ingestion of hot gases into 
the inlet. The other mechanism of hot gas ingestion, which is 
not as severe, is due to the recirculating flow caused by the 
interaction of the outward flowing hot gases with the 
headwind. The outward flowing hot gases are forced to 
stagnate by the freestream flow and are deflected towards the 
engine inlet where they are reingested. 

Figure 3 shows the velocity distribution in an x-y plane passing 
through the center of the fore and aft jets for the four cases (for 


the sake of clarity, different vector scales are used in different 
segments of the domain). The effect of headwind to jet velocity 
ratio (Uo«/Uj) on the flow field is clearly evident In Case 1 and 
2 (U«yUj=0.03), the effect of the lift jets extends far upstream 
(35 jet diameters for Case 1 and 46 jet diameters for Case 2), 
while for Case 3 and 4 (Uoo/Uj=0.09) the stagnation region is 
much closer to the jets. For the low velocity ratio, hot gas from 
the lift jets flow under the fuselage towards the inlet, where part 
of it is ingested directly by the strong suction of the inlet, while 
the rest is carried by its streamwise momentum further 
upstream until it is forced to stagnate by the oncoming flow and 
recirculate back towards the inlet (ground vortex). The most 
striking difference between high and low U«AJj is the complete 
absence of the recirculating ground vortex in Case 3 and 4. 
This is due to the high forward momentum of the headwind 
which prevents the hot gas from the jets from moving 
upstream. The effect of different h/d ratios is not as strong. For 
h/d of 2 (Case 2 and Case 4), a smaller area is available 
between the fuselage and the ground and consequently the 
radial momentum of the hot gas is larger than for h/d of 4. This 
enhances the spread of the jet, both in the lateral and 
streamwise direction. A comparison of Case 3 and Case 4 (just 
upstream of the fore jet and under the fuselage), suggests that 
for h/d of 2 the forward flow is stronger than for h/d of 4. This 
is also the case for Case 1 and 2, in which the recirculating 
zone for h/d of 2 extends about 1 1 jet diameters further 
upstream than for h/d of 4. 

Figure 4 shows temperature contours corresponding to the 
velocity vector fields shown in Figure 3. In all cases, the 
temperature fields upstream of the lift jets are very similar to 
those implied by the velocity field. For Cases 1 and 2 
(U«/Uj=0.03), the hot gases from the jets are transported much 
further upstream of the inlet Although not clearly evident in 
Figure 4, the upstream transport and the subsequent 
recirculation of hot gases (although cooler) back to the inlet can 
be directly correlated with the velocity vector field. This effect 
can be clearly seen in Figure 5, which summarizes Figures 3 
and 4. In Figure 5(a) and 5(b), the cooler recirculated gases 
(positive streamwise velocity) can be clearly distinguished from 
the hotter gases (negative streamwise velocity). This relation 
between the velocity field and the temperature field is also 
evident in Case 4 (U«/Uj=0.09, h/d = 2) between the inlet and 
the fore jet where hot gases are transported upstream along the 
bottom of the fuselage. However, for Case 3 (U«yUj=0.09, h/d 
= 4), the temperature field extends further upstream than 
indicated by the velocity field; this indicates either a diffusive 
transport of heat or convective transport of heat from the sides. 
Another interesting feature is the difference in temperature 
between the fore half of the fountain and the aft half. In all 
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cases the fore half is cooler than the aft half, which must be due 
to the its greater interaction with the cooler freestream fluid. 
This effect was also observed in the simulations of [16]. 
Maclean et al. [20] compared the upstream distribution of 
temperature (along the jet center and along the axis of 
symmetry) obtained from their experiments with the present 
calculations for Case 3. In the region near the jets, the 
agreement between the two is good, but further upstream the 
experiments show a higher concentration of jet fluid. 

Figure 6 shows the velocity distribution in the ground plane for 
Case 1 (Uoo/Uj=0.03, h/d = 4) and Case 3 (Uoo/Uj=0.09, h/d = 
4). In the vicinity of the jets, the outward flowing radial wall 
jets are clearly visible. These wall jets interact with each other 
to form the streamwise and spanwise fountain flows. We can 
also see the effect of the headwind on the outward movement of 
the hot jet fluid; for the low velocity ratio the momentum of the 
jets carry the hot gases much further than with the higher 
headwind. This effect is more clearly seen in Figures 7 and 8, 
which show line and color contours of the temperature field in 
the ground plane. For the low headwind case, both the forward 
and lateral spread of the hot gases is much more than with the 
higher headwind. In addition, the height of the fuselage above 
the ground (h/d ratio) also has a effect on the forward and 
lateral spreading of the jets. As mentioned earlier, for h/d =2 the 
radial momentum of the jet is higher and this can be seen in the 
greater lateral spread of the jet fluid before it is pushed 
downstream by the headwind. Another interesting feature is the 
sharp temperature gradients observed between the jets in the 
fountain upwash region. This is due to the vertical movement 
of two jet streams at different temperatures with little transverse 
movement of fluid across this boundary. In fact, the locus of 
points across which these sharp gradients are seen in the lateral 
direction is the stagnation line which divides fluid from the fore 
and aft jets until finally they mix and smooth out the 
differences. These sharp gradients are also present in the 
simulations of [16] and in the experiments of [20], although the 
experimentally observed gradients arc not as sharp. 

Table 2 compares the upstream extent of the recirculating hot 
gas zone with the experiments of [20]. In the present study, the 
temperature distribution is used as a means to determine the 
upstream spread of the hot gases (Figures 4 and 5). It is seen 
that the calculated values agree well with the experimental 
values for Cases 1 , 3 and 4, while for Case 2 the calculated 
value of x/d is much higher. Also, the trend indicated by the 
experiments (increased upstream penetration (x/d) with 
increased h/d ratio) is the opposite of what the numerical 
calculations indicate. This in effect would indicate that the jet 
fluid mixes more rapidly with the free stream flow and loses its 


Table 2. Extent of Recirculating Hot Gas Zone 
Upstream of Foreward Lift Jets. 


CASE 

Present 

Calculations 

Experiments 
Ref. [20] 

1 

35 

31 

2 

46 

22 

3 

10 

12 

4 

12.5 

14 


identity earlier than the calculations predict. This discrepency 
could be due to the inability of the k-E turbulence model to 
accurately capture the full extent of turbulence production 
caused by the stagnating jets. Recent numerical studies of the 
fountain region of twin jet impingement [28] indicate that the 
turbulent kinetic energy is severely under predicted, particularly 
in the region close to the ground plane. This in pan also 
explains the sharp temperature gradients seen in the upwash 
region of the fountain flow. The authors suspect that for Case 2 
this might be a critical factor in determining the dilution and 
spread of the jet fluid. 

Figures 9-12 show velocity and temperature distributions in 
three z-y planes for each of the cases tested. Pan (a) is the z-y 
plane passing through the fountain upwash between the jets, 
part (b) is the z-y plane between the fore jet and the inlet and 
part (c) is the plane at the inlet. In all four cases, the spanwise 
fountain flow which is responsible for the movement of hot 
gases upstream from the fore jets is clearly visible in (b). As 
the hot gases from the jets move upstream under the fuselage 
they are considerably diluted due to mixing with the cooler 
headwind. As this gas approaches the inlet, pan of it is directly 
ingested into the inlet (c). It is interesting to see that for the low 
velocity ratio, ingestion occurs from the side, top and bottom of 
the fuselage while for the higher headwind most of the 
ingestion occurs from the bottom and side of the fuselage. 
Another interesting feature is the movement of hot gases; for 
the higher headwind, the hot gases show a lateral outward and 
upward movement while for the lower headwind the flow is 
laterally inward, i.e., cooler fluid is entrained into the fountain 
region. This phenomenon is also reflected in the ground plane 
temperatures (Figures 7 and 8), where the lateral spread of the 
jet fluid between the fore jet and the inlet is more with the 
higher headwind. This is due to the higher dynamic pressure of 
the oncoming free stream, which forces the hot gases out from 
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Table 3. Non-dimensional Mass Averaged Inlet 

Plane Temperatures. (T ayg . T„)/(Tj - T„) 


CASE 

Present 

Calculations 

Caculations 
Ref. [16] 

1 

0.083 

0.13 

2 

0.097 

0.11 

3 

0.0 

0.08 

4 

0.017 

0.03 


under the fuselage. Table 3 shows the mass averaged non- 
dimensional inlet plane temperatures calculated for the four 
cases and compares them to those obtained by the numerical 
simulations of [16]. There is some discrepancy between the 
two calculations. It is not possible to identify the reasons for 
this discrepancy, but the difference in the resolution of the flow 
field between the two methods could be one of the major 
contributing factors. The results of the present calculation 
follow the trends indicated by the velocity and the temperature 
profiles. The maximum average inlet temperature is for Case 2 
(Uoo/Uj=0.03, h/d = 2) while Case 3 does not seem to ingest 
any hot gases. 

4. Summary and Conclusions 

In the present paper, a multigrid calculation procedure for three- 
dimensional flows is used to study the hot gas ingestion by 
STOVL aircraft. Calculations with a simulated fuselage and 
twin exhaust jets have been made for two ground positions and 
two head-wind ratios. The global features of the flowfield like 
the ground plane temperature distributions and the formation of 
the ground vortex for the four cases are according to 
expectations and agree well with the recent numerical study of 
ref. 16. A major contribution of the present study is the 
demonstration that the multigrid algorithm can significantly 
reduce the computational effort required to solve the governing 
equations compared to a conventional single grid procedure. 
The reduction in computational effort permits more parametric 
studies to be conducted. 

In the present study, the Reynolds-averaged flow equations 
have been solved in conjunction with the k-e turbulence model. 
The flow is considered to be steady in time, although, it has 
been observed [27] that the region of jet impingment is highly 


unsteady and the turbulence structures are very complex. 
Therefore, the accuracy of the current calculations is limited by 
the assumptions made in deriving the governing equations of 
the turbulence model. A study performed on the twin jet 
impingement flow in isolation [28] showed that the k-e model 
is able to predict the velocity field reasonably well, but the 
turbulence kinetic energies are severely under-predicted. 
Although it is possible to incorporate more complex turbulence 
models into the current algorithm, at this time it is not certain 
that they improve the predictions over those of the k-e model. 
Future extension of this study will be in the direction of 
representing the aircraft geometry more realistically through the 
use of curvilinear grids. 
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Figure 2: Typical convergence history for multigrid calculation procedure. 100 x 44 x 48 
cells with 2 coarse levels. 
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Figure 3: Velocity vectors in x-y plane passing through the jet centerline ( z= 0)- Plot 
domain extends from -46 <x/d^ 14 and from the ground plane to y/d - 2U. 
Vector scales vary in different regions of the flow field. (a) Case 1 (h/d — , 
Uoc/Uj = 0.03); (b) Case 2 (h/d = 2, Uoo/Uj = 0.03); (c) Case 3 (h/d = 4, Uoo/Uj 
= 0.09); (d) Case 4 (h/d = 2, Uoo/Uj = 0.09). 
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Figure 4: 


Temperature contours in x-y plane passing through the jet centerline (z=0). 
Contour levels: - - - 300-475 K; 475-650 K; . 650-825 K; 


825-1000 K. (a) Case 1 (h/d = 4, Uco/Uj = 0.03); (b) Case 2 (h/d=2, 

UWUi = 0.03); (c) Case 3 (h/d = 4, Uoo/Uj = 0.09); (d) Case 4 (h/d=2, U«/Uj 

= 0.09). 
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(d) 


Figure 5. Superimposed velocity vectors and temperature color contours in x-y plane 

through the jet centerline (Figures 3 and 4). Plot domain extends from -46 < x/d 
^14 and from the ground plane to y/d=20.(a) Case 1; (b) Case 2; (c) Case 3- 
(d) Case 4. 
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Figure 6: Velocity vectors in x-z ground plane. Plot domain extends from -10 < x/d < 14 
and from symmetry plane to z/d=12.(a) Case 1 (h/d=4, Uoo/Uj = 0.03); (b) Case 
3 (h/d=4, Uoo/Uj = 0.09). 
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Figure 7: Temperature contours in x-z ground plane. Contour levels: — 300-475 K; 475-650 K; . _ 

650-825 K; 825-1000 K. (a) Case 1 (h/d=4, U^yUj = 0.03); (b) Case 2 (h/d=2, Uoo/Uj=0.03); (c) 
Case 3 (h/d=4, U«/Uj = 0.09); (d) Case 4 (h/d=2, U«/Uj=0.09). 









Figure 8: Temperature color contours in x-z ground plane (Figure7). (a) Case 1; (b) Case 
2; (c) Case 3; (d) Case 4. 
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Figure 9: 


Velocity vector and temperature contours in different z-y planes for Case 1 (h/d 
= 4 Uoo/Ui = 0.03). Plot domain extends from symmetry plane to z/d - y ana 
from ground plane to y/d = 5. (a) x/d = 3 (between fore an a t jets in e _ 
fountain upwash region); (b) x/d = -5 (between inlet and fore j ), ( ) 

( at inlet). 
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Figure 11: Velocity vector and temperature contours in different z-y planes for Case 3 (h/d 
= 4, IWUj = 0.09). Plot domain extends from symmetry plane to z/d = 9 and 
from ground plane to y/d = 5. (a) x/d = 3 (between fore and aft jets in the 
fountain upwash region); (b) x/d = -5 (between inlet and fore jet); (c) x/d = -10 
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Figure 12: Velocity vector and temperature contours in different z-y planes for Case 4 (h/d 
= 2, Uoo/Uj = 0.09). Plot domain extends from symmetry plane to z/d = 9 and 
from ground plane to y/d = 5. (a) x/d = 3 (between fore and aft jets in the 
fountain upwash region); (b) x/d = -5 (between inlet and fore jet); (c) x/d = -10 
(at inlet). 
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Ingestion of hot exhaust gases by the engines of Short Take-off and Vertical Landing (STOVL) aircraft has 
been an important research problem for several years. The hot gas environment around STOVL aircraft is 
three-dimensional and turbulent. In this study, the Navier-Stokes equations governing the hot gas ingestion 
flowfield are solved by an efficient finite-difference calculation procedure. The complete geometry, including 
the head wind and the fuselage, is simulated. Four demonstration calculations with variations in the height of 
the fuselage and the head wind velocity are presented. It is shown that the calculation procedure efficiently 
provides a solution to the governing equations, and produces realistic descriptions of the flowfield and temper- 
ature field. 


Introduction 

A NUMBER of interesting fluid dynamic effects have been 
identified when the lift jets of Short Take-off and Ver- 
tical Landing (STOVL) aircraft impinge on the ground sur- 
face. 1-3 First, as the lift jets expand, they entrain the sur- 
rounding fluid, causing a negative pressure underneath the 
fuselage and a loss in lift. As the jets impinge on the ground 
and spread radially outwards, the wall jets further entrain 
external fluid and increase the loss in lift. In a multiple- jet 
configuration, these wall jets collide with each other and turn 
upwards to form an upwash fountain. This fountain flow has 
two main effects on the STOVL aircraft dynamics. First, an 
increase in lift force is caused when the fountain impinges on 
the aircraft fuselage. The recovery in lift is a positive effect 
of the upwash flow. However, this impinging fluid can also 
flow along the fuselage surface and eventually make its way 
into the engine inlets. Because the temperature of the fountain 
flow is much hotter than the ambient air, its ingestion by the 
engine can reduce the power and cause thermal stresses in 
the components. In addition to the fountain flow, another 
mechanism for hot gas ingestion results from the interaction 
of the forward-moving wall jet with the headwind. When the 
headwind and the wall jet collide, a stagnation region is formed, 
and the wall jet is turned into a ground vortex 2 This ground 
vortex can subsequently be re-ingested by the engine, result- 
ing in a further loss in power. However, the hot gas ingestion 
process depends on where and how the inlets to the engine 
are located with respect to the approach flow. For strong 
headwinds, the ground vortex can be pushed behind the en- 
gine inlets, such that no hot gases are ingested. It can be 
recognized that the overall flowfield governing these fluid 
dynamic interactions is strongly three-dimensional, as well as 
turbulent. 

In the last three decades, a number of studies addressing 
these issues has been conducted (see Ref. 3 for good review). 
A majority of these studies were experimental in nature, al- 
though some analytical studies of the individual processes, 
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such as single jet impingement, have also been conducted. 
These studies have identified the fundamental processes and, 
to a certain extent, quantified the effect of the parameters on 
hot gas ingestion and fountain flow. However, detailed mea- 
surements of the velocity and temperature fields in the foun- 
tain region and in the ground vortex in realistic operating 
conditions are not currently available. Recently, Saripalli 4 
conducted some experiments in a water tunnel in which Laser- 
Doppler velorimeter and laser imaging techniques were used 
to study the fountain region of the lift jets. These studies 
involved isothermal flows, and were limited to the fountain 
formation. 

With recent advances in the development of powerful dig- 
ital computers, it has now become possible to numerically 
solve the equations that govern the transport of mass, mo- 
mentum, and enthalpy. However, because the flow is tur- 
bulent, such studies must necessarily make assumptions on 
the macroscopic behavior of turbulence. Despite the inherent 
limitations of such an approach in representing all the effects 
of turbulence with precision, it has become a useful engi- 
neering tool. Numerical computation of the complete STOVL 
aircraft flowfield requires the solution of nonlinear partial- 
differential equations that govern the transport of mass, mo- 
mentum, and heat in three-dimensional space with boundary 
conditions that describe the lift jets, aircraft fuselage, and the 
head wind. Because of the complexity of such solutions and 
limited computer capabilities, until recently most numerical 
studies were limited to two dimensions, in which simplifica- 
tions in the spanwise direction were made. Kotansky and 
Bower 5 were probably the first investigators to apply a Navier- 
Stokes analysis to the study of planar jet impingement of 
relevance to Vertical Take-Off and Landing (VTOL) aircraft. 
In their analysis, a one-equation model for turbulence kinetic 
energy was integrated with equations describing the stream 
function and transport of vorticity. A numerical solution of 
an impinging jet was obtained and compared with experi- 
mental data. This study was followed by Agarwal and Bower, 6 
who improved the turbulence model by considering an ad- 
ditional equation for the turbulence dissipation. Because the 
improved model explicitly calculated the local-length scale of 
turbulence, the flowfields predicted by their model displayed 
better agreement with experimental data. 

Although the hot gas ingestion process results primarily 
from the impingement of the lift jets, the eventual ingestion 
is dictated by the complete flow dynamics, including the for- 
mation of the fountain upwash and its interaction with the 
headwind. For stronger headwind speeds, it is possible to push 
the ground vortex downstream of the engine inlet and, thus. 
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completely avoid its contribution to hot gas ingestion. Also, 
the collision of the wall jets and the fountain formation can 
be affected by the headwind and the geometry of the fuselage 
under surface. In addition, placement of Lift Increasing De- 
vices (LIDs ) 2 and other obstructions can divert the fountain 
flow from the engine inlets, thereby reducing the severity of 
the hot gas ingestion. In order to better understand all the 
flow processes, a complete Navier-Stokes analysis, including 
the far-field condition of a prescribed headwind, is necessary. 
Although such an analysis is computationally very intensive, 
the benefits are significant for evaluating the effects of geo- 
metric and flow parameters. Recently, an important contri- 
bution has been made in this direction by VanOverbeke and 
Holdeman . 7 In this study, the complete fuselage, headwind, 
and multiple impinging jets were numerically simulated, and 
the temperature fields close to the engine inlet face were 
studied. By varying the head wind and the impingement height, 
they were able to alter the ingestion pattern and quantify their 
effects. A finite-difference computer program* was used, and 
the flow domain was discretized into a relatively large number 
of finite volumes. 

A primary drawback of numerical computations of three- 
dimensional elliptic fluid flows is that they can require sub- 
stantial amounts of computer time to obtain a converged so- 
lution. These large computational costs often discourage sys- 
tematic and thorough investigations of the effects of various 
parameters. However, with research into faster converging 
numerical algorithms, it is possible to substantially reduce 
these computer times. In this paper, we demonstrate the ap- 
plicability of one such algorithm , 9 which is based on the con- 
cept of using multiple levels of grids 10 to obtain faster con- 
vergence. In this study, the geometrical configuration studied 
by VanOverbeke and Holdeman 7 is considered, and a set of 
four calculations have been performed. To compare our re- 
sults with those of VanOverbeke and Holdeman , 7 we have 
simulated exactly the same configurations, albeit with a finer 
grid, and compared the computer times and the flowfields. 
The four calculations in this study differ in the height of the 
jets from the ground and the ratio of the head wind to the 
jet velocity. 

The following sections briefly describe the solution pro- 
cedure and present the results of the calculations. In Part 1 
of this study, a separate experimental investigation is pre- 
sented by McLean et al., n in which the concentration of the 
exhaust jets is quantified through a laser imaging technique. 
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Governing Equations and Solution Algorithm 

The flowfield and temperature field around the STOVL 
aircraft are governed by the following set of elliptic partial- 
differential equations. In the present analysis, the flow is con- 
sidered to be steady in the time-averaged sense, and the ef- 
fects of turbulence are captured by a turbulence model . 12 The 
governing equations can be written as follows: 
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where u, v y w, are the velocities in the x, y, and z directions, 
respectively; P is the pressure, R is the gas constant, and G 
is the production of turbulent kinetic energy. In all these 
equations, the diffusive coefficients contain the effects of lam- 
inar as well as turbulent diffusion. Standard values of the 
constants in the turbulence model 12 are used. 

These equations are solved on a rectangular domain that 
envelopes the aircraft fuselage and the exhaust jets. Because 
the current computer program is limited to the use of rectan- 
gular grids, the fuselage of the aircraft was modeled to be of 
rectangular shape. However, efforts are currently under way 
to extend the program to handle grids of arbitrary inclinations 
(boundary-fitted grids), which will then permit a more realistic 
representation of the aircraft shape and angles of attack. The 
boundary conditions and solution domain considered in cur- 
rent calculations will be described in a later section. 

The above equations are finite-differenced by a hybrid 
scheme 9 combining first- and second-order differencing for 
the total convective and diffusive operator. At high-cell Peclet 
numbers, the finite-differencing scheme becomes first-order 
accurate to preserve stability of the solution procedure. The 
multigrid solution technique 10 differs from the single grid tech- 
nique in the manner the finite -difference equations are solved 
to the required accuracy. It is well-known that, for elliptic 
equations, single grid techniques converge poorly when the 
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finite-difference grid contains a large number of mesh points. 
This results from the low-frequency errors that are slow to 
converge on a given fine grid. In the multigrid concept, a 
series of fine and coarse grids are used, and the solution is 
switched between the coarse and fine grids such that errors 
of all frequencies converge at the same rate. This novel con- 
cept was originally proposed for model elliptic equations, but 
recently it has become popular in solving practical fluid flows. 

The principle behind the multigrid procedure may be ex- 
plained as follows. Consider an elliptic equation that is dis- 
cretized by finite-difference or finite element methods. Gen- 
erally, an iterative procedure that is designed to obtain a 
solution to the discrete equations will converge rapidly in the 
first few iterations, but subsequently the convergence dete- 
riorates. It can be shown through formal analyses that the 
cause for this poor convergence is the presence of low-fre- 
quency errors that are not resolved well on any given grid. 
However, these errors can converge faster if the grid is coarse; 
therefore, in the multigrid procedure, these errors are inter- 
polated to a coarser grid and solved. The results of this coarse 
grid solution are then put back on the finer grid solution 
through extrapolation. In the case of a calculation with many 
mesh points, several layers of coarse meshes can be used to 
obtain rapid convergence. However, the coarsest grid is often 
dictated by the constraints of geometry and boundary con- 
ditions. 

In the present solution scheme, we have combined the mul- 
tigrid technique with a coupled solution scheme for the mo- 
mentum and continuity equations. The velocities and pres- 
sures are obtained by a block-implicit solution of the momentum 
and continuity equations. Because this procedure is explained 
in previous works by one of the authors, 9 it is not elaborated 
on here. However, in order to simulate the fuselage and the 
lift jets, the procedure has been extended to handle internal 
obstacles and mass and momentum injections. These features 
have been incorporated in a manner that maintains consist- 
ency between the coarse and fine grids — a necessary feature 
for the algorithm to eventually converge to the correct so- 
lution. It must be pointed out that, although the coarse grids 
are used, the final solution is obtained on the finest grid in 
the system. Thus the intermediate grids are used only for 
accelerating the convergence rate, and do not influence the 
final converged solution. 

The present computer program was initially validated in 
model laminar flow problems prior to being used for more 
complex flows. Recently, the computer program was applied 
to calculate the interaction of a transverse jet with a cross- 
stream flow. 13 Also, the program was validated by applying 
it to studying the twin-jet upwash flow in isolation and 
comparing 14 the resulting solutions with the experimental data 
of Saripalli. 4 These two studies, as well as a calculation of a 
single jet impingement (unpublished), provide confidence in 
the use of the present computer program to the complete 
STOVL configuration. 

The next section describes the computational domain and 
a typical convergence history of the multi grid solution pro- 
cedure. Simulated flowfields and temperature fields are pre- 
sented, along with a discussion on the mechanism of hot gas 
ingestion for a four- jet configuration. These results compare 
qualitatively with the flow imaging data of McLean, Sullivan, 
and Murthy, 11 as well as with calculations presented by Van 
Overbeke and Holdeman. 7 

Results 

Computational Details and Convergence Histories 

In this study, four test cases have been considered, which 
are similar to those studied in Ref. 7. The important flow and 
geometrical parameters are outlined in Table 1 . The four cases 
differ in the height of the fuselage from the ground plane (h/ 
d, where d is the diameter of the lift jets) and the headwind 
to jet velocity ratio (UJU t ). The solution domain is a three- 


Tabk 1 Flow and geometry parameters 


Case 

UJU, 

h/d 

Ma, 

-» 

o 

X 

oc 

Re, x 10- s 

1 

0.03 

4.0 

0.47 

3.3 

3.3 

2 

0.03 

2.0 

0.47 

3.3 

3.3 

3 

0.09 

4.0 

0.47 

9.9 

3.3 

4 

0.09 

2.0 

0.47 

9.9 

3.3 



a) 



Fig. 1 Computational domain and grid distribution: a) side view; b) 
front view (not to same scale as a». 

dimensional wind tunnel with overall dimensions of 136d in 
the streamwise direction (x), 27 or 29d in the cross-stream 
direction (y), and 40d in the spanwise direction (z). These 
large dimensions are necessary to accurately resolve the flow- 
field around the jets and the fuselage without any interference 
from the simulated boundaries. Figure 1 shows the compu- 
tational grid used and the flow geometry. The geometry sim- 
ulates four jets (due to symmetry, only half the domain is 
considered) issuing from the engine with the engine inlet lo- 
cated 10 jet diameters upstream of the fore jet axis. The 
distance between the fore and aft jet axis is 6 diameters; the 
spanwise distance between the jets is 3.2 diameters. The fu- 
selage extends throughout the calculation domain, and is 2.5 
diameters in height and 2.2 diameters in width. The inlet, 
which is 67 jet diameters from the entrance of the test section, 
has a height of 2.5 diameters and a width of 1.4 diameters. 
It is adjacent to the fuselage in the spanwise direction and 
extends 9.5 jet diameters up to the fore jet, after which it 
becomes part of the fuselage. The calculation domain is di- 
vided into 100 computational cells in the x direction, 44 in 
the y , and 48 in the z direction, giving a total of about 21 1 ,000 
cells. The particular grid distribution was chosen to provide 
good resolution in the near field of the jets and, at the same 
time, keep the grid expansion ratio as small as possible. The 
grid distribution was also influenced by the necessity to ac- 
commodate internal obstacles and baffles, such that grid lines 
bounded the surfaces of these structures. 

At the inlet of the test section, uniform velocity and tem- 
perature conditions are prescribed for the incoming head- 
wind; the turbulent kinetic energy was assumed to be 1% of 
the mean kinetic energy. At the exit, outflow boundary con- 
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ditions are prescribed. In they direction, wall boundary con- 
ditions are prescribed. The bottom wall is assumed to have 
zero heat loss, whereas the top wall is assumed to be at a 
uniform temperature. In the span wise direction, symmetry 
boundary conditions are applied on one side, whereas an 
adiabatic wall is prescribed on the other side. Within the 
calculation domain, additional conditions are prescribed for 
the engine inlet and the jets. The inlet acts as a mass sink. 



Fig. 2 Typical convergence history for multigrid calculation proce- 
dure, 100 x 44 x 48 cells with two coarse levels. 


the jets as mass sources. In prescribing these conditions, it is 
important to conserve mass. The jet velocity profile is as- 
sumed to be uniform, with the turbulent kinetic energy equal 
to 1% of the mean kinetic energy. 

In the calculation procedure, the prescribed fine grid was 
coarsened to two additional levels. The calculation was started 
on the coarsest grid from prescribed initial conditions, but, 
on the finer grids, the starting solutions were obtained by 
successively extrapolating converged flowfields from the coarser 
grids. This Full Multigrid Cycle (FMG) provided realistic 
starting values for the calculations on the fine grids. Figure 2 
shows a typical convergence history for the calculation of case 
1. The convergence is monitored by the normalized sum of 
the absolute mass and momentum residuals over the whole 
solution domain. It is seen that the solution converges to the 
required tolerance of 1% error in approximately 150 itera- 
tions. In the previous study of VanOverbeke and Holdeman, 7 
which used a single grid procedure, the number of iterations 
were approximately 2000 for a 134,000 grid. This represents 
a significant reduction in computational effort. The present 
calculation required approximately 35 min of CPU time on a 
CRAY-2; typically, this translates to a few hours of interactive 
real clock time. 


Calculated Velocity and Temperature Fields 

Hot gas ingestion is a serious problem for STOVL aircraft 
in ground proximity. Two important factors governing hot gas 



Fig. 3 Superimposed velocity vectors and temperature gray scale contours in x-y plane through the jet centerline ( z = 0). Plot domain extends 
from —46 ^ x/d ^ 14 and from the ground plane to y/d = 20: a) case 1; b) case 2; c) case 3; d) case 4; headwind from left to right. 
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ingestion are the proximity of the lift jets to the ground (h/d 
ratio) and the ratio of the freestream velocity to the jet ve- 
locity {UJU } ).- Two other important factors, not considered 
in the present study are the location of the inlet with respect 
to the life jets and the use of shields and baffles to deflect 
the flow of hot gases away from the engine inlet. When the 
lift jets impinge on the ground, they stagnate and form radially 
flowing wall jets. When two opposing wall jets meet, they 
manifest themselves into a fountain flow. The flowfield caused 
by a rectangular four-jet configuration is governed by a strong 
fountain flow between the fore and aft jets (streamwise foun- 
tain) and also between the symmetrical spanwise jets (span- 
wise fountain). The streamwise fountain upwash, when ob- 
structed by the fuselage, is forced to flow laterally outwards 
and up, whereas the spanwise fountain upwash tends to flow 
under and along the length of the fuselage in the streamwise 
direction. Both these mechanisms contribute to the direct 
ingestion of hot gases into the inlet. The other mechanism of 
hot gas ingestion, which is not as severe, is due to the recir- 
culating flow caused by the interaction of the outward flowing 
hot gases with the headwind. The outward flowing hot gases 
are forced to stagnate by the freestream flow, and are de- 
flected towards the engine inlet where they are reingested. 

Figure 3 shows velocity vectors superimposed on gray scale 
temperature contours in the x-y plane passing through the 
center of the fore and aft jets for the four cases. (For the sake 
of clarity, different vector scales are used in different segments 
of the domain.) The effect of headwind to jet velocity ratio 
(UJU f ) on the flowfield is clearly evident. In cases 1 and 2 
(UJU f = 0.03), the effect of the lift jets extends far upstream 
(35 jet diameters for case 1 and 46 jet diameters for case 2), 
whereas for cases 3 and 4 (UJU f = 0.09), the stagnation 
region is much closer to the jets. For the low-velocity ratio, 
hot gas from the lift jets flows under the fuselage towards the 
inlet, where part of it is ingested directly by the strong suction 
of the inlet, while the rest is carried by its streamwise mo- 
mentum further upstream until it is forced to stagnate by the 
oncoming flow and recirculate back towards the inlet (ground 
vortex). The most striking difference between high and low 
UJU t is the complete absence of the recirculating ground 
vortex in cases 3 and 4. This is due to the high forward mo- 
mentum of the headwind, which prevents the hot gas from 
the jets from moving upstream. The effect of different h/d 
ratios is not as strong. For h/d of 2 (cases 2 and 4), a smaller 
area is available between the fuselage and the ground and, 
consequently, the radial momentum of the hot gas is larger 
than for h/d of 4. This enhances the spread of the jet, both 
in the lateral and streamwise direction. A comparison of cases 
3 and 4 (just upstream of the fore jet and under the fuselage) 
suggests that for h/d of 2 the forward flow is stronger than 
for h/d of 4. This is also the case for cases 1 and 2, in which 
the recirculating zone for h/d of 2 extends about 11 jet di- 
ameters further upstream than for h/d of 4. 

In all cases, the temperature fields upstream of the lift jets 
are very similar to those implied by the velocity field. For 
cases 1 and 2 ( UJU ) = 0.03), the hot gases from the jets are 
transported much further upstream of the inlet. The upstream 
transport and the subsequent recirculation of hot gases (al- 
though cooler) back to the inlet can be directly correlated 
with the velocity vector field. In Figs. 3a and 3b, the cooler 
recirculated gases (positive streamwise velocity) can be clearly 
distinguished from the hotter gases (negative streamwise ve- 
locity). This relation between the velocity field and the tem- 
perature field is also evident in case 4 (UJUj = 0.09, h/d = 
2) between the inlet and the fore jet, where hot gases are 
transported upstream along the bottom of the fuselage. How- 
ever, for case 3 ( UJU i = 0.09, h/d = 4), the temperature 
field extends further upstream than indicated by the velocity 
field. This indicates either a diffusive transport of heat or 
convective transport of heat from the sides. Another inter- 
esting feature is the difference in temperature between the 
fore half of the fountain and the aft half. In all cases, the fore 


half is cooler than the aft half, which must be due to its greater 
interaction with the cooler freestream fluid. This effect was 
also observed in the simulations in Ref. 7. Maclean et al. n 
compared the upstream distribution of temperature (along 
the jet center and along the axis of symmetry) obtained from 
their experiments with the present calculations for case 3. In 
the region near the jets, the agreement between the two is 
good, but further upstream the experiments show a higher 
concentration of jet fluid. 

Figure 4 shows gray scale contours of the temperature field 
in the ground plane. For the low headwind case, there is much 
more forward and lateral spread of the hot gases than with 
the higher headwind. In addition, the height of the fuselage 
above the ground (h/d ratio) also has an effect on the forward 
and lateral spreading of the jets. As mentioned earlier, for 
h/d = 2, the radial momentum of the jet is higher, as can be 
seen in the greater lateral spread of the jet fluid before it is 
pushed downstream by the headwind. Another interesting 
feature is the sharp temperature gradients observed between 
the jets in the fountain upwash region. This is due to the 
vertical movement of two jet streams at different tempera- 
tures, with little transverse movement of fluid across this 
boundary. In fact, the locus of points across which these sharp 
gradients are seen in the lateral direction is the stagnation 
line that divides fluid from the fore and aft jets, until finally 



Fig. 4 Temperature contours in x-z ground plane: a) case 1; b) case 
2; c) case 3; d) case 4; headwind from left to right. 


Table 2 Extent of recirculating hot gas zone upstream of 
foreword lift jets 


Case 

Present 

Calculation 

Experiments, 
Ref. 11 

1 

35 

31 

2 

46 

22 

3 

10 

12 

4 

12.5 

14 
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c) 


Fig. 5 Velocity vector and temperature contours in different z-y planes for case 1. Plot domain extends from symmetry plane to zld - 9, and 
from ground plane to y! d — 5: a) x/d = 3; b) x/d = — 5; c) x/d — — 10. 


they mix and smooth out the differences. These sharp gra- 
dients are also present in the simulations in Ref. 7 and in the 
experiments in Ref. 11 , although the experimentally observed 
gradients are not as sharp. 

Table 2 compares the upstream extent of the recirculating 
hot gas zone with the experiments in Ref. 11. In the present 
study, the temperature distribution is used as a means to 
determine the upstream spread of the hot gases (Fig. 3). It is 
seen that the calculated values agree well with the experi- 
mental values for cases 1, 3, and 4, whereas for case 2 the 
calculated value of x/d is much higher. Also, the trend indi- 


cated by the experiments [increased upstream penetration 
(x/d) with increased h/d ratio] is the opposite of what the 
numerical calculations indicate. This, in effect, would indicate 
that the jet fluid mixes more rapidly with the freestream flow 
and loses its identity earlier than the calculations predict. This 
discrepancy could be due to the inability of the k-s turbulence 
model to accurately capture the full extent of turbulence pro- 
duction caused by the stagnating jets. Recent numerical stud- 
ies of the fountain region of twin jet impingement 14 indicate 
that the turbulent kinetic energy is severely underpredicted, 
particularly in the region close to the ground plane. This, in 


47 









D. K. TAFTI AND S. P. VANKA 


J. AIRCRAFT 




Temperature Levels 

300 - 425 K 

... 425 • 550 K 

550 • 675 K 

~ 675 - *00 K 



b) 



Fig. 6 Velocity vector and temperature contours in different z-y planes for case 3. Plot domain extends from symmetry plane to z/d = 9, and 
from ground plane toy/d = 5: a) x/d = 3; b) x/d = -5; c) x/d = - 10. 


part, also explains the sharp temperature gradients seen in 
the upwash region of the fountain flow. The authors suspect 
that, for case 2, this might be a critical factor in determining 
the dilution and spread of the jet fluid. More comparisons 
between the present simulations and experiments can be found 
in Ref. 11. 

Figures 5 and 6 show velocity and temperature distributions 
in three z-y planes for cases 1 and 3. Part a) is the z-y plane 
passing through the fountain upwash between the jets, part 

b) , the z-y plane between the fore jet and the inlet, and part 

c) , the plane at the inlet. In both cases, the spanwise fountain 


flow responsible for the movement of hot gases upstream from 
the fore jets is clearly visible in b). As the hot gases from the 
jets move upstream under the fuselage, they are considerably 
diluted, due to mixing with the cooler headwind. As this gas 
approaches the inlet, part of it is directly ingested into the 
inlet c). It is interesting to see that, for the low-velocity ratio, 
ingestion occurs from the side, top, and bottom of the fuse- 
lage, whereas, for the higher headwind, most of the ingestion 
occurs from the bottom and side of the fuselage. Another 
interesting feature is the movement of hot gases: for the higher 
headwind, the hot gases show a lateral outward and upward 
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Table 3 Nondimensional mass averaged inlet plane temperatures 
(F. vg - TMTj - TJ 


Case 

Present 

Calculation 

Calculation, 
Ref. 7 

1 

0.083 

0.13 

2 

0.097 

0.11 

3 

0.0 

0.08 

4 

0.017 

0.03 


movement whereas for the lower headwind the flow is laterally 
inward, i.e. , cooler fluid is entrained into the fountain region. 
This phenomenon is also reflected in the ground plane tem- 
peratures (Fig. 4), where the lateral spread of the jet fluid 
between the fore jet and the inlet is greater with the higher 
headwind. This is due to the higher dynamic pressure of the 
oncoming freestream, which forces the hot gases out from 
under the fuselage. Table 3 shows the mass averaged non- 
dimensional inlet plane temperatures calculated for the four 
cases, and compares them to those obtained by the numerical 
simulations in Ref. 7. There is some discrepancy between the 
two calculations. It is not possible to identify the reasons for 
this discrepancy, but the difference in the resolution of the 
flowfield between the two methods could be one of the major 
contributing factors. The results of the present calculation 
follow the trends indicated by the velocity and the tempera- 
ture profiles. The maximum average inlet temperature is for 
case 2 {UJUj = 0.03, h/d = 2); case 3 does not seem to 
ingest any hot gases. 

Summary and Conclusions 

In the present paper, a multigrid calculation procedure for 
three-dimensional flows is used to study the hot gas ingestion 
by STOVL aircraft. Calculations with a simulated fuselage 
and twin exhaust jets have been made for two ground posi- 
tions and two head- wind ratios. The global features of the 
flowfield (e.g., the ground plane temperature distributions 
and the formation of the ground vortex) for the four cases 
are according to expectations and agree well with the recent 
numerical study of Ref. 7. A major contribution of the present 
study is the demonstration that the multigrid algorithm can 
significantly reduce the computational effort required to solve 
the governing equations, compared to a conventional single 
grid procedure. The reduction in computational effort permits 
more parametric studies to be conducted. 

In the present study, the Reynolds-averaged flow equations 
have been solved in conjunction with the k-e turbulence model. 
The flow is considered to be steady in time, although, it has 
been observed 11 that the region of jet impingement is highly 
unsteady and the turbulence structures are very complex. 
Therefore, the accuracy of the current calculations is limited 
by the assumptions made in deriving the governing equations 
of the turbulence model. A study performed on the twin jet 
impingement flow in isolation 14 showed that the k-e model is 
able to predict the velocity field reasonably well, but the 
turbulence kinetic energies are severely underpredicted. Al- 
though it is possible to incorporate more complex turbulence 


models into the current algorithm, it is not certain at this time 
that they improve the predictions over those of the k-e model. 
Future extension of this study will be in the direction of rep- 
resenting the aircraft geometry more realistically through the 
use of curvilinear grids. 
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ABSTRACT 


Hot gas ingestion problems for STOVL 
(Short Take-Off, Vertical Landing) aircraft are 
typically approached with empirical methods and 
experience. In this study, the hot gas environment 
around a STOVL aircraft was modeled as multiple 
jets in crossflow with inlet suction. The flow field 
was calculated with a Navier-Stokes, Reynolds- 
averaged, turbulent, 3D CFD code using a multi- 
grid technique. A simple model of a STOVL aircraft 
with four choked jets at 1000 K was studied at 
various heights, headwind speeds, and thrust splay 
angles in a modest parametric study. Scientific 
visualization of the computed flow field shows a pair 
of vortices in front of the inlet. This and other 
qualitative aspects of the flow field agree well with 
experimental data. 

NOMENCLATURE 


D ( = characteristic length of jet nozzles, 0.0366 
m (1.44 in) 

H = distance from ground to bottom of the 
aircraft (aircraft altitude or height) 

U_ = headwind velocity 

V, = jet nozzle exit velocity, 633 m/s (2080 ft/sec) 
8 = thrust splay angle measured from the 

downward vertical toward symmetry plane 
k = turbulent kinetic energy 
e = turbulent energy dissipation 
x = axial Cartesian coordinate, zero at upstream 
boundary 

x-y = vertical plane aligned in axial direction 
x-z = horizontal plane 

y = vertical Cartesian coordinate, zero at ground 
plane 

y-z = vertical plane aligned in spanwise direction 


z = spanwise Cartesian coordinate, zero at 
aircraft centerline plane 

INTRODUCTION 


Hot gas ingestion can cause significant 
problems for a STOVL (Short Take-Off, Vertical 
Landing) aircraft including reduced thrust and 
compressor stalls. These problems involve many 
hazards for the pitots including very hard landings. 
During the design of a 
ingestion problems are 
empirical methods and 
power of today’s supercomputers and workstations, 
numerical methods employing efficient algorithms 
are becoming a viable engineering tool for analysis 
and design. In a previous endeavor, VanOverbeke 
& Holdeman 34 proved the feasibility of CFD analysis 
for hot gas ingestion. This study is a follow-on 
effort exploring the practicality of using an efficient 
numerical method for the problem of hot gas 
ingestion. 


STOVL aircraft, hot gas 
typically approached with 
experience . Given the 


FLOW FIELD DESCRIPTION 


Ingestion of hot gases generates problems 
in two ways: an average temperature rise results in 
a loss of engine thrust, and a temperature distortion 
may cause the engine to stall. Engine exhaust 
gases may be ingested by far-field and/or near-field 
mechanisms. A schematic of these mechanisms is 
shown in Fig. 1. 

The far-field mechanism results from the 
exhaust gases impinging on the ground and forming 
radial wall jets which flow forward, separate, and 
mix with the headwind. Near-field ingestion occurs 
with multiple jet configurations. Wall jets flowing out 
from the lift jets meet and create an upflow or 
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Hot gas ingestion mechanisms 

Figure 1 


fountain. This fountain flow can impinge on the 
aircraft’s underside, flow along the fuselage to the 
engine inlets, and be ingested. The gases ingested 
by this near-field mechanism tend to be hotter, 
giving greater temperature distortion than those 
ingested by the far-field mechanism. 

As stated earlier, the hot gas environment 
around a STOVL aircraft was modeled as multiple 
jets in crossflow with inlet suction. Mass sources 
represent the nozzle exits, and a mass sink at the 
end of the inlet provides the suction. The mass 
injected by the nozzles balances the mass removed 
by the inlet suction. This configuration derives from 
the previous study by VanOverbeke & Holdeman 34 . 
To meet the requirements of the CFD code, the 
aircraft model is placed in a confined flow, i.e., a 
’wind tunnel’. Also, the aircraft model has no angle 
of attack due to the use of a cartesian grid based 
flow solver. 

The STOVL aircraft model (see Fig. 2) is 
composed of rectangular solids for the fuselage and 
engine. For computational simplicity, the nose and 



the tail of the aircraft reach to infinity, and the 
model lacks wings. Baffles on the sides of the 
fuselage comprise the walls of the inlet. The 
nozzles are square in cross-section and are flush 
with the bottom of the aircraft. The square cross- 
section of the jets and the rectangular aircraft body 
result from the use of the cartesian grid. 

The four choked nozzles inject air at 1000 K 


(1340°F) straight down into the flowfield with a 
velocity of 633 m/s (2080 ft/sec). Each lift jet 
issues from the nozzle exit in a uniform flow. The 
headwind is also a uniform flow, but at a 
temperature of 300 K (81 °F). This approximates an 
aircraft landing with a forward speed or an aircraft 
facing into a wind which lacks a boundary layer. 
In the baseline case, the headwind (U_) flows at 3% 
of the jet velocity (V), or about 19 m/s (37 kts), and 
the distance from the ground to the bottom of the 
aircraft model (H) is four times the characteristic 
length of the nozzles (4 D). The parametric studies 
include: various altitudes (H = 2 to 32 D) at a 
constant headwind (U„ = 0.03 V,); various headwind 
speeds (U_ = 0.01 to 0.09 V,) for a constant aircraft 
altitude (H = 4 D,); and various thrust splay angles 
(5 = 0° to 45°) for a constant height (H = 4 Dj) 
and headwind speed (U„ = 0.03 V,). 

The physical dimensions of the aircraft 
model are given in Table I. Note that the forward 
and aft nozzles have the same side-to-side 
separation, i.e., they are in-line, not offset. 


Table I 

STOVL aircraft model dimensions 


Fuselage: 


width (nose) = 

2.25 D, 

width (tail) = 

5.05 D, 

height * 

2.5 D, 

length = 

oo 

Inlet: 


width 

1.4 D, 

height = 

2.5 D, 

length = 

9.5 D, 

Jets separation: 


(center-to-center) 


side-to-side = 

3.25 D : 

fore & aft = 

6.0 D, 


NUMERICAL DESCRIPTION 
Calculation domain 

The grid geometry used for the baseline 
case is shown in Fig. 3. Exhibited are the 
centerline plane, the ground plane, and a vertical 
spanwise plane at the end of the domain as well as 
the aircraft model. The grid shows the high density 
of the calculation nodes in the region of the jets. 
For all calculations, symmetry assumptions allowed 
calculating only half of the physical domain. 

The other boundary conditions for the 
calculation domain include an inflow simulating the 
headwind for the domain face in front of the aircraft 
model and an outflow condition for the domain face 
behind the model. All flow properties are defined 
for the inflow condition. The outflow condition, in 
contrast, merely assumes the properties of the 
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Figure 3 



Grid for headwind speed variation 


axially nearest cells. The 
top, bottom, and remaining 
side of the domain are no- 
slip, stationary walls as are 
the aircraft surfaces. All the 
walls assume adiabatic 
conditions for the energy 
equation. The symmetry 
plane also has a symmetry 
condition for the energy 
equation. 

For the height 
variation, the grid contains 
211,200 cells arranged as 
follows: 100 cells in the x 
direction, 44 cells in the y 
direction, and 48 cells in the 
z direction. The physical 
dimensions of the baseline 
grid are about 135 Q long, 

29 Dj high, and 40 Dj wide. 

The aircraft to ground 
distance was varied by 
elongating the cells 
underneath the aircraft. This 
facilitated the comparative 
analysis. 

The headwind speed 
variation used a slightly 
modified grid; 100 cells in x, 

44 cells in y, and 60 cells in 
z yielding 264,000 total cells. 

(See Fig. 4.) This grid has 
a greater length (177 Dj) and 
greater width (59 Dj) than the 
baseline grid. Also, the 
distance in front of the 
forward pair of jets is greater (152 Dj versus 76 D, 
for the baseline grid) to accommodate the long 
region of hot gas in front of the inlet in the U_ = 
0.01 Vj case. No grid modification was needed to 
vary the headwind speed. 

Flow solver 

the flow field in this domain was calculated 
with a Navier-Stokes, Reynolds-averaged CFD code. 
References 5-7 describe this steady-state CFD code 
(and its techniques) for the three-dimensional 
analysis of turbulent elliptic flows in a Cartesian 
coordinate system. 

The CFD code (CART3D) solves the time- 
averaged Navier-Stokes or Reynolds equations. 
The k-e turbulence model provides closure. The 
governing equations include continuity, x-, y-, and z- 
momenta, energy, turbulent kinetic energy, and 
turbulent energy dissipation. These equations are 
solved using a block-implicit multigrid algorithm 
developed by Vanka. 

CART3D uses a hybrid differencing scheme 
on a staggered grid. This means that the code 


Figure 4 

uses central differencing or upwind differencing 
depending upon the cell’s Reynolds number. Also, 
the scalar properties (density, pressure, etc.) are 
calculated at the cell volume centers while the 
velocities are solved at the centers of the cell faces. 

The multi-grid technique speeds 
convergence by solving the equations on sequential 
grids of different cell densities. The flow is 
initialized on the coarsest grid which gets refined by 
the multi-gridding. Dividing each cell on a grid into 
eight equal cells refines the grid for the next grid 
level. A V-cycle of sweeps on the various grid 
levels is performed until the solution converges on 
the finest grid. This technique speeds convergence 
by dampening out errors with the various levels of 
grid refinement. 

To determine convergence, the residuals are 
non-dimensionalized by an appropriate number, and 
then the maximum of all the residuals is compared 
to the tolerance criterion. The tolerance criterion 
used by this study is 1% for the finest grid. All test 
cases used the third grid level for the finest grid. 
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FLOW FIELD FEATURES 

A short study of the 
features in the baseline case 
will help bring out the 
differences caused by varying 
the aircraft altitude. 

Fig. 5 displays the 
temperature contours in an x- 
z plane near the ground. 

These contours show the 
locations of the forward 
vortex pair and the two 
ground vortices generated by 
the interaction of the jets and 
the crossflow. The axis for 
the forward vortex pair is 
perpendicular to the plane 
shown in Fig. 5 while the 
axes for the ground vortices 
are parallel to the plane. 

The forward vortex pair is 
smeared by the steady-state 
calculations but still agrees 
well with the time-averaged 
experimental data 8 ' 10 . 

The particle traces in 
Fig. 6 reveal ingestion of 
exhaust gases. The particle 
traces from the jet region 
show the forward vortex and 
ingestion into the inlet. The 
particle traces starting at the 
inflow boundary show the 
headwind’s deflection around 
the forward vortex. 

Temperature contours along 
the ground plane are also 
shown for clarity of position. 

Fig. 7 shows a three- 
dimensional temperature 
contour for the baseline 
case. This isotherm is for 
325 K (125°F), a reasonable 
upper limit on the 
temperature of the air 
reaching the engine. With 
the ambient flow at 300 K 
(81 °F), this represents a 
temperature rise of more 
than 25 K (45°F) for the fluid 
inside the isotherm. The 
inlet is almost completely 
obscured by the contour. 

Clearly, the engine is 
exposed to a considerable 
amount of hot gas from the 
engine exhaust. The bubble in front of the inlet 7 is due to wall effects which have no consequence 

reveals the location of the forward vortex. Note that on the hot gas ingestion, 

the clipping of the isotherm in the right side of Fig. 



H - 4 D„ U. - 0.03 V., T_ - 300 K, T, - 1000 K 

Figure 7 



Select particle traces 
H - 4 0., U_ - 0.03 V, 

Figure 6 



x-z plane temperature contours 
H - 4 D., U_ - 0.03 V, 

Figure 5 
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AIRCRAFT ALTITUDE 
PARAMETER VARIATION 

The aircraft altitude varied from 2 D, to 32 
Dj. The cases actually computed over this range (H 
= 2, 3, 4, 5, 6, 8. 12, 16, 24, and 32 D.) were 
carefully chosen to capture the changes in the flow 
field features. Fig. 8 shows temperature contours 
for select aircraft altitudes (2 D,, 4 Q, 8 Dj, 12 Dj, 
and 32 D,) in two planes: the ground plane and a 
vertical plane passing through the jets near the jet 
centers. These cases show the changes in the flow 
field affecting hot gas ingestion over the range of 
variation. In each case, the fuselage is mostly 
hidden by the vertical plane of temperature 
contours. The major effects of aircraft altitude can 
be seen in this figure: the forward vortex changes 
in character, the amount of hot gas ingested is 
reduced, and the ground vortices decrease in size. 

In Fig. 9, the temperatures at the cells in 
front of the mass sink or 'engine face' are plotted 
against the aircraft altitude. The temperatures 
shown are the minimum, the maximum, and a 
weighted average based on the cell volumes. The 
spread of the minimum temperature and the 
maximum temperature shows the temperature 
distortion at the engine face. This and the average 
temperature rise above ambient is plotted in Fig. 
10. At low altitudes, the distortion is obviously 



Temperature contours in select planes 
(U. - 0.03 V^ 

Figure 8 



0 5 to 19 20 25 30 35 

Aircraft Altitude * Dj 


Plot of engine face temperatures 
vs. aircraft altitude 

Figure 9 

extreme; and it quickly diminishes with increasing 
altitude. The average temperature shows a similar 
behavior. The non-monotonic behavior at H = 3 D, 
in all of these curves appears to be physical. Other 
calculations on different grids exhibit the same 
trends. 



0 5 10 15 20 25 3C 35 


Aircraft Altitude - Dj 


Distortion and average rise 
vs. aircraft altitude 

Figure 10 

HEADWIND SPEED PARAMETER VARIATION 

The headwind speed was varied from 1 % to 
9% of the jet velocity (U. = 0.01 to 0.09 V) in 1% 
increments for a constant altitude (H = 4 Dj). Fig. 
1 1 displays temperature contours in the same two 
planes as in Fig. 8: the ground plane and a vertical 
plane passing through the jets near the jet centers. 
Selected headwind speeds (1%, 3%, 5%, 7%, and 
9%) show the changes in the flow field for the 
varying headwind speed. Again, the fuselage is 
mostly hidden by the vertical plane of temperature 
contours. In Fig. 1 1 the calculated domain includes 
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Temperature contours in select planes 


(H - 4 D,) 

Figure 11 


unrealistically representing a windy vertical landing, 
but they might be relevant for a low speed landing. 
For choked jets, the 9% headwind represents about 
a 110 knot headwind which would either be a 
hurricane or a slower than normal landing speed for 
a conventional fighter aircraft. 


H.4D, 
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Distortion and average rise 
vs. headwind speed 
Figure 13 


a much larger region in front of the inlet in 
comparison with Fig. 8. This is due to a very long 
region of exhaust gas which extends in front of the 
inlet along the ground in the U_ = 0.01 V. case. 

The effect of headwind speed on tne engine 
face temperatures can be seen in Fig. 12. Note 
that the minimum temperature declines rather 
steadily with increasing headwind speed. The 
temperature distortion at the engine face varies 
weakly with headwind speed at the low speeds and 
is greatest for the 7% case as shown in Fig. 13. 
One should note that these high velocity ratios are 


H - 4 D, 
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Plot of engine face temperatures 
vs. headwind speed 

Figure 12 


THRUST SPLAY ANGLE PARAMETER VARIATION 

A technique used to help control hot gas 
ingestion is to splay the jets. By angling the lift 
jets, the relative strengths of the fountain, upwash, 
and vortices are changed thus changing the flow 
structures affecting ingestion. In this study, the 
splay angle (8) of the thrust is measured from the 
downward vertical inward to the centerline plane of 
the aircraft model. To vary the thrust splay angle, 
the component velocities on the jets changed to 
provide the required angle while keeping the speed 
of the jet constant. Thus the direction of the lift jets 
changed while the geometry of the aircraft model 
did not change. 

Splaying all lets 

For the first variation of the thrust splay 
angle, all four jets were splayed the same amount. 
The splay angle varied from 0* to 45' (8 = 0' to 
45’) in 5' increments for a constant height (H = 4 
Dj) and constant headwind speed (U_ = 0.03 V,). 
Fig. 14 shows temperature contours in two planes 
similar to those displayed in Fig. 8 & Fig. 11: the 
ground plane and a vertical plane that is now on 
the inner side of the lift jets (instead of near the 
center). Selected splay angles (O', 5', 10’, 20', 
25’, and 30') show the changes in the flow 
structures due to varying the angle of all the jets. 
The most noticeable change occurs in the length of 
the hot gas region in front of the inlet. A less 
noticeable change is the increase of hot gas directly 
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Temp. contours * splay all jets 
(H - 4 D j( U. - 0,03 V,) 

Figure 14 


in front of the inlet for the 5" case. This actually 
causes an increase in hot gas ingestion over the 0° 
case, instead of reducing it. 

This and other effects of thrust splay angle 
for all the jets can be seen in Fig. 15 and Fig. 16. 
In Fig. 15, the average temperature first rises and 
then drops until a splay angle of about 20" where 
it starts rising again. A local maximum exists at the 
25* point before the average temperature flattens 



Plot of engine face temperatures 
vs. splay angle on all jets 
Figure 15 


out at its lowest value. Both the minimum and 
maximum temperatures follow the same behavior. 
The distortion and average temperature rise 
displayed in Fig. 16 show the same patterns. Note 
that for this configuration of model and altitude, the 
jets will converge at the ground for a thrust splay 
angle of 22*. 
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Distortion and average rise 
vs splay angle on all jets 

Figure 16 



Splaying forward lets only 

In the second variation, only the forward jets 
varied in thrust splay angle. The rear lift jets 
maintained a 0’ splay angle. Again, the thrust 
splay angle varied from 0* to 45" (S = 0" to 45") in 
5’ increments at the same altitude (H = 4 D,) and 
headwind speed (U_ = 0.03 V.) just as in the ail jets 
splayed variation. Fig. 17 snows the temperature 
contours in the same two planes as in Fig. 14: the 
ground plane and a vertical plane on the inner side 
of the lift jets. Again, the selected splay angles (O’, 
5", 10*, 20*. 25’, and 30") show the changes in the 
flow structures due to varying the angle of the 
forward jets alone. The flow field changes are 
basically the same as in the cases with all the jets 
splayed. 

The ingestion effects of thrust splay angle 
for the forward jets can be seen in Fig. 18. The 
average temperature first rises and then drops until 
a splay angle of about 20* where it rises very 
slightly. A local maximum exists at 25" splay angle 
before the average temperature drops and flattens 
out at its lowest value. Both the minimum and 
maximum temperatures follow a similar behavior. 
The distortion and average temperature rise 
displayed in Fig. 19 show the same patterns, just 
as when all the jets are splayed. Overall, splaying 
the forward jets alone gives the same effects on hot 
gas ingestion as splaying all the jets. 
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Forward jets splayed 
(H - 4 Dj, U.. - 0.03 V,) 

Figure 17 


CONCLUSIONS/REMARKS 

In summary, the average of the engine face 
temperatures decreases with increasing height but 
is relatively unaffected by headwind speed. Engine 
face temperature distortion also decreases with 
increasing height, but unexpectedly increases with 
headwind speed until the forward vortex is behind 
the inlet face. The headwind speed variation 
reveals that (for a constant height) the engine face 
temperature is dominated by near field ingestion 
effects. 

As for the thrust splay angle variation, 
splaying the jets inward (for a constant height and 
headwind speed) first causes a rise and then a 
rapid decrease in hot gas ingestion, although a 
local maximum exists when the jets converge at the 
ground plane. Also, splaying the forward jets alone 
instead of all the jets gives almost the same 
reduction in hot gas ingestion without as large a 
thrust penalty. 

A comparison of the H = 4 Dj, U = 0.03 V, 
cases in the height parameter variation and the 
headwind speed parameter variation show the 
effects of the calculation domain on the flow field. 
The vertical walls of the narrower calculation 
domain definitely affect the ground vortices, but no 
differences exist in the temperatures reaching the 



vs. splay angle on forward jets 

Figure 18 


H « 4 D,, U„ « 0.03 V, 



Distortion and average rise 
vs. splay angle on forward jets 

Figure 10 


engine. Essentially, the tunnel walls in the 
calculation domain used for the height parameter 
variation are sufficiently far from the aircraft model 
so as to not affect the desired quantities (engine 
face temperatures). If the overall flow field is of 
primary interest, then the tunnel walls would have 
to be farther from the aircraft model. 

This study did not address the importance 
of the aircraft geometry (fuselage, wings, tails, etc.) 
in relation to the flow field. Only one aircraft model 
was used, and it was quite simplistic. 

The last conclusion from this study concerns 
the practicality of using an efficient CFD code for 
parameter variation studies. The turn-around time 
on a Cray-2 supercomputer and state-of-the-art 
workstations allows quick parameter changes. 
Typically, a Cray-2 supercomputer solved the flow 
field in about an hour with a turn-around time of a 
day. A dedicated IBM RS-6000 workstation can 
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solve the flow field in about 6 hours and can 
actually give shorter turn-around than the shared 
supercomputer. 
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with empirical methods and experience, in this study, the hot gas environment around a STOVL aircraft was 
modeled as multiple jets in crossflow with inlet suction. The flowfield was calculated with a Navier-Stokes, 
Reynolds-averaged, turbulent, three-dimensional CFD code using a multigrid technique. A simple model of a 
STOVL aircraft with four choked jets at 1000 K was studied at various heights, head wind speeds, and thrust 
splay angles in a modest parametric study. Scientific visualization of the computed flowfield shows a pair of 
vortices in front of the inlet. 


Nomenclature 

D - characteristic length of jet nozzles, 0.0366 m 
(1.44 in.) 

H = distance from ground to bottom of the aircraft 
(aircraft altitude or height) 
k = turbulent kinetic energy 
t/» = head wind velocity 
y. = jet nozzle exit velocity, 633 m/s (2080 ft/s) 
x = axial Cartesian coordinate, zero at upstream 
boundary 

x-y = vertical plane aligned in axial direction 
x-z = horizontal plane 

y = vertical Cartesian coordinate, zero at ground plane 
y-z = vertical plane aligned in spanwise direction 
z — spanwise Cartesian coordinate, zero at aircraft 
centerline plane 

5 = thrust splay angle measured from the downward 

vertical toward symmetry plane 
e = turbulent energy dissipation 

Introduction 

H OT gas ingestion can cause significant problems for a 
STOVL (short take-off, vertical landing) aircraft, in- 
cluding reduced thrust and compressor stalls. These problems 
involve many hazards for the pilots, including very hard land- 
ings. During the design of a STOVL aircraft, hot gas ingestion 
problems are typically approached with empirical methods 
and experience. 1 - 2 Given the power of today’s supercomputers 
and workstations, numerical methods employing efficient al- 
gorithms are becoming a viable engineering tool for analysis 
and design. In a previous endeavor, VanOverbeke and 
Holdeman 3 - 4 demonstrated the feasibility of CFD analysis for 
hot gas ingestion. This study is a follow-on effort exploring 
the practicality of using an efficient numerical method for the 
problem of hot gas ingestion. 
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Flowfield Description 

Ingestion of hot gases generates problems in two ways: 1) 
an average temperature rise results in a loss of engine thrust, 
and 2) a temperature distortion may cause the engine to stall. 
Engine exhaust gases may be ingested by far-field and/or near- 
field mechanisms. A schematic of these mechanisms is shown 
in Fig. 1. 

The far-field mechanism results from the exhaust gases im- 
pinging on the ground and forming radial wall jets which flow 
forward, separate, and mix with the head wind. Near-field 
ingestion occurs with multiple jet configurations. Wall jets 
flowing out from the lift jets meet and create an upflow or 
fountain. This fountain flow can impinge on the aircraft’s 
underside, flow along the fuselage to the engine inlets, and 
be ingested. The gases ingested by this near-field mechanism 
tend to be hotter, giving greater temperature distortion than 
those ingested by the far-field mechanism. 

As stated earlier, the hot gas environment around a STOVL 
aircraft was modeled as multiple jets in crossflow with inlet 
suction. Mass sources represent the nozzle exits, and a mass 



Fig. 1 Hot gas ingestion mechanisms. 



(HEADWIND) AND TEMPERATURE SPECIFIED 
Fig. 2 STOVL Aircraft model. 
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sink at the end of the inlet provides the suction. The mass 
injected by the nozzles balances the mass removed by the 
inlet suction. This configuration is derived from the previous 
study by VanOverbeke and Holdeman. 34 To meet the re- 
quirements of the CFD code, the aircraft model is placed in 
a confined flow, i.e., a “wind tunnel.” (It should be noted 
that changes in the width of the confinement did not signif- 
icantly affect the flowfield in the vicinity of the jets and inlet 
for the same inflow and jet conditions. Thus, the wind-tunnel 
walls had minimal effect on the desired results.) Also, the 
aircraft model has no angle of attack due to the use of a 
Cartesian grid-based flow solver. 

The STOVL aircraft model (see Fig. 2) is composed of 
rectangular solids for the fuselage and engine. For compu- 
tational simplicity, the nose and the tail of the aircraft reach 
to infinity, and the model lacks wings. Baffles on the sides of 
the fuselage comprise the walls of the inlet. The nozzles are 
square in cross section and are flush with the bottom of the 
aircraft. The square cross section of the jets and the rectan- 
gular aircraft body result from the use of the Cartesian grid. 

The four choked nozzles inject air at 1000 K (1340°F) straight 
down into the flowfield with a velocity of 633 m/s (2080 ft/s). 
Each lift jet issues from the nozzle exit in a uniform flow. 
The head wind is also a uniform flow, but at a temperature 
of 300 K (81°F). This approximates an aircraft landing with 


Table 1 STOVL Aircraft 
model dimensions 


Fuselage 

Width (nose) 

2.25D ; 

Width (tail) 

5.05D / 

Height 

2.5D, 

Length 

30 

Inlet 

Width 

\AD, 

Height 

2.5 D, 

Length 

9.5 D, 

Jets separation 
Center-to-center 

Side-to-side 

3.25 D } 

Fore and aft 

6.0D, 


a forward speed, or an aircraft facing into a wind which lacks 
a boundary layer. In the baseline case, the head wind (U x ) 
flows at 3% of the jet velocity (V,), or about 19 m/s (37 kt), 
and the distance from the ground to the bottom of the aircraft 
model (H) is four times the characteristic length of the nozzles 
(4 Dj). The parametric studies include: various altitudes (H 
= 2-32D ; ) at a constant head wind (£/* = 0.03 V y ); various 
head wind speeds ( U * = 0.01-0. 09 V y ) for a constant aircraft 
altitude ( H = 4 D)\ and various thrust splay angles (8 = 0- 
45 deg) for a constant height ( H = 4 D f ) and head wind speed 

(u x = 0.03 vy. 

The physical dimensions of the aircraft model are given in 
Table 1. Note that the forward and aft nozzles have the same 
side-to-side separation, i.e., they are in-line, not offset. 

Numerical Description 

Calculation Domain 

The grid geometry used for the baseline case is shown in 
Fig. 3. Exhibited are the centerline plane, the ground plane, 
and a vertical spanwise plane at the end of the domain as well 
as the aircraft model. The grid shows the high density of the 
calculation nodes in the region of the jets. For all calculations, 
symmetry assumptions allowed calculating only half of the 
physical domain. 

The other boundary conditions for the calculation domain 
include an inflow simulating the head wind for the domain 
face in front of the aircraft model and an outflow condition 
for the domain face behind the model. All flow properties 
are defined for the inflow condition. The outflow condition, 
in contrast, merely assumes the properties of the axially near- 
est cells. The top, bottom, and remaining side of the domain 
are no-slip, stationary walls, as are the aircraft surfaces. All 
the walls assume adiabatic conditions for the energy equation. 
The symmetry plane also has a symmetry condition for the 
energy equation. 

For the height variation, the grid contains 211,200 cells 
arranged as follows: 100 cells in the x (longitudinal) direction, 
44 cells in the y (vertical) direction, and 48 cells in the z 
(horizontal) direction. The physical dimensions of the base- 
line grid are about 135Z) y long, 29 D } high, and 40D y wide. The 


VERTICAL SPANWISE PLANE 


CENTERLINE PLANE 



GROUND PLANE 


CENTERLINE PLANE 


\ NOZZLE 
INLET 

FUSELAGE 

Fig. 3 Grid for baseline case. 


VERTICAL SPANWISE PLANE 



NOZZLE 


GROUND PLANE 


Fig. 4 Grid for head wind speed variation. 
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aircraft to ground distance was varied by elongating the 
cells underneath the aircraft. This facilitated the comparative 
analysis. 

The head wind speed variation used a slightly modified grid; 
100 cells in x, 44 cells in y, and 60 cells in z, yielding 264,000 
total cells (see Fig. 4). This grid has a greater length (177D,) 
and greater width (59 D.) than the baseline grid. Also, the 
distance in front of the forward pair of jets is greater (152D / 
vs 76 D i for the baseline grid) to accommodate the long region 
of hot gas in front of the inlet in the £/„ = 0.01 V) case. (In 
each of the grids, the inflow boundary condition is far enough 
from the stagnation point of the hot ground flow to prevent 
interference with the fluid mechanics. Also, the tunnel wall 
boundaries are sufficiently far from the aircraft model to min- 
imize effects on the flowfield.) No grid modification was needed 
to vary the head wind speed. 

Flow Solver 

The flowfield in this domain was calculated with a Navier- 
Stokes, Reynolds-averaged CFD code. References 5-7 de- 
scribe this steady-state CFD code (and its techniques) for the 
three-dimensional analysis of turbulent elliptic flows in a 
Cartesian coordinate system. 

The CFD code (CART3D) solves the time-averaged Na- 
vier-Stokes or Reynolds equations. The k-e turbulence model 
provides closure. The governing equations include continuity, 
x-, y-, and z-momenta, energy, turbulent kinetic energy, and 
turbulent energy dissipation. These equations are solved using 
a block-implicit multigrid algorithm developed by Vanka. 5 

CART3D uses a hybrid-differencing scheme on a staggered 
grid. This means that the code uses central differencing or 
upwind differencing, depending upon the cell’s Reynolds 
number. Also, the scalar properties (density, pressure, etc.) 
are calculated at the cell volume centers, while the velocities 
are solved at the centers of the cell faces. 

The multigrid technique speeds convergence by solving the 
equations on sequential grids of different cell densities. The 
flow is initialized on the coarsest grid which gets refined by 
the multigridding. Dividing each cell on a grid into eight equal 


cells refines the grid for the next grid level. A V cycle of 
sweeps on the various grid levels is performed until the so- 
lution converges on the finest grid. This technique speeds 
convergence by dampening out errors with the various levels 
of grid refinement. 

To determine convergence, the residuals are nondimen- 
sionalized by an appropriate number, and then the maximum 
of all the residuals is compared to the tolerance criterion. The 
tolerance criterion used by this study is 1% for the finest grid. 
All test cases used the third grid level for the finest grid. 

Flowfield Features 

A short study of the features in the baseline case will help 
bring out the differences caused by varying the aircraft alti- 
tude, the head wind speed, and the thrust splay angle. 

Figure 5 displays the temperature contours in an x-z plane 
near the ground. These contours show the locations of the 
forward vortex pair and the two ground vortices generated 
by the interaction of the jets and the crossflow. The axis for 
the forward vortex pair is perpendicular to the plane shown 
in Fig. 5, while the axes for the ground vortices are parallel 
to the plane. The forward vortex pair is smeared by the steady- 
state calculations, but still agrees well with the time-averaged 
experimental data. 8-10 

The particle traces in Fig. 6 reveal ingestion of exhaust 
gases. The particle traces from the jet region show the forward 
vortex and ingestion into the inlet. The particle traces starting 
at the inflow boundary show the head wind’s deflection around 
the forward vortex. Temperature contours along the ground 
plane are also shown for clarity of position. 

Figure 7 shows a three-dimensional temperature contour 
for the baseline case. This isotherm is for 325 K (125°F), a 
reasonable upper limit on the temperature of the air reaching 
the engine. With the ambient flow at 300 K (81°F), this rep- 
resents a temperature rise of more than 25 K (45°F) for the 
fluid inside the isotherm. The inlet is almost completely ob- 
scured by the contour. Clearly, the engine is exposed to a 
considerable amount of hot gas from the engine exhaust. The 



Fig. 5 x-z plane temperature contours H = 4 D j% U ^ = 0.03V . 
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GROUND PLANE TEMPERATURE CONTOURS 

Fig. 7 Three-dimensional isotherm (T = 325 K) H = 4 D Jt V * = 0.03V,, T. = 300 K, Tj = 1000 K. 



Fig. 8 Temperature contours in select planes (V x = 0.03V,). 

bubble in front of the inlet Veveals the location of the forward 
vortex. Note that the clipping of the isotherm in the right side 
of Fig. 7 is due to wall effects which have no consequence on 
the hot gas ingestion. 



Aircraft Altitude - Dj 


Fig. 9 Plot of engine face temperatures vs aircraft altitude 




Fig. 10 Distortion and average rise vs aircraft altitude. 


Aircraft Altitude Parameter Variation 

The aircraft altitude varied from 2 to 32 D r The cases ac- 
tually computed over this range ( H = 2, 3, 4, 5, 6, 8, 12, 16, 
24, and 32D,) were carefully chosen to capture the changes 
in the flowfield features. Figure 8 shows temperature contours 
for select aircraft altitudes (2, 4, 8, 12, and 32D y ) in two 
planes: 1) the ground plane and 2) a vertical plane passing 
through the jets near the jet centers. These cases show the 
changes in the flowfield affecting hot gas ingestion over the 
range of variation. In each case, the fuselage is mostly hidden 
by the vertical plane of temperature contours. The major 
effects of aircraft altitude can be seen in this figure: the for- 
ward vortex changes in character, the amount of hot gas in- 
gested is reduced, and the ground vortices decrease in size. 

In Fig. 9, the temperatures at the cells in front of the mass 
sink or “engine face” are plotted against the aircraft altitude. 
The temperatures shown are the minimum, the maximum, 
and a weighted average based on the cell volumes. The spread 
of the minimum and maximum temperature shows the tem- 
perature distortion at the engine face. This, and the average 
temperature rise above ambient, is plotted in Fig. 10. At low 
altitudes, the distortion is obviously extreme, and it quickly 
diminishes with increasing altitude. The average temperature 
shows a similar behavior. The nonmonotonic behavior at H 
= 3 D } in all of these curves appears to be physical. At this 
time, the physical mechanism creating the local increase at H 


= 3 D i is unknown. In tests for grid dependencies, other cal- 
culations on different grids exhibited the same trends, in- 
cluding the nonmonotonicity. 

Head Wind Speed Parameter Variation 

The head wind speed was varied from 1 to 9% of the jet 
velocity - 0.01-0. 09 V,) in 1% increments for a constant 
altitude (H — 4D / ). Figure 11 displays temperature contours 
in the same two planes as in Fig. 8; the ground plane and a 
vertical plane passing through the jets near the jet centers. 
Selected head wind speeds (1, 3, 5, 7, and 9%) show the 
changes in the flowfield for the varying head wind speed. 
Again, the fuselage is mostly hidden by the vertical plane of 
temperature contours. In Fig. 11 the calculated domain in- 
cludes a much larger region in front of the inlet in comparison 
with Fig. 8. This is due to a very long region of exhaust gas 
which extends in front of the inlet along the ground in the U * 
= 0.01 V, case. 

The effect of head wind speed on the engine face temper- 
atures can be seen in Fig. 12. Note that the minimum tem- 
perature declines rather steadily with increasing head wind 
speed. The general flatness of the average temperature would 
indicate that near-field ingestion dominates over most of the 
speed range. The temperature distortion at the engine face 
varies weakly with head wind speed at the low speeds, and is 
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Fig. 14 Temperature contours in select planes — splaying all jets 
(H = 4 D Jy U x = 0.03V,). 
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Fig. 12 Plot of engine face temperatures vs head wind speed. 



Headwind Speed • % of Vj 

Fig. 13 Distortion and average rise vs head wind speed. 

greatest for the 7% case as shown in Fig. 13. One should note 
that these high-velocity ratios are unrealistically representing 
a windy vertical landing, but they might be relevant for a low- 
speed runway landing. For choked jets, the 9% head wind 
represents about a 110-kt head wind which would either be 
a hurricane, or a slower than normal landing speed for a 
conventional fighter aircraft. 

Thrust Splay Angle Parameter Variation 

A technique used to help control hot gas ingestion is to 
splay the jets. By angling the lift jets, the relative strengths 


— = 3% — * 4 

V| q 



K 250 1 i 1 1 1 1 1 i 1 i i I 

-5 0 5 1 0 1 5 20 25 30 35 40 45 50 

Splay Angle on Ail Jets - deg. 

Fig. 15 Plot of engine face temperatures vs splay angle on all jets. 

of the fountain, upwash, and vortices are changed, therefore 
changing the flow structures affecting ingestion. In this study, 
the splay angle (5) of the thrust is measured from the down- 
ward vertical inward to the centerline plane of the aircraft 
model. To vary the thrust splay angle, the component veloc- 
ities on the jets changed to provide the required angle while 
keeping the speed of the jet constant. Thus, the direction of 
the lift jets changed while the geometry of the aircraft model 
did not. 

Splaying All Jets 

For the first variation of the thrust splay angle, all four jets 
were splayed the same amount. The splay angle varied from 
0 to 45 deg (5 = 0-45 deg) in 5-deg increments for a constant 
height ( H = 4 D f ) and constant head wind speed = 0.03 V,). 
Figure 14 shows temperature contours in two planes similar 
to those displayed in Figs. 8 and 11; the ground plane and a 
vertical plane that is now on the inner side of the lift jets 
(instead of near the center). Selected splay angles (0, 5, 10, 
20, 25, and 30 deg) show the changes in the flow structures 
due to varying the angle of all the jets. The most noticeable 
change occurs in the length of the hot gas region in front of 
the inlet. A less noticeable change is the increase of hot gas 
directly in front of the inlet for the 5-deg case. This actually 
causes an increase in hot gas ingestion over the 0-deg case, 
instead of reducing it. 

This and other effects of thrust splay angle for all the jets 
can be seen in Figs. 15 and 16. In Fig. 15, the average tem- 
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perature first rises and then drops until a splay angle of about 
20 deg, where it starts rising again. A local maximum exists 
at the 25-deg point before the average temperature flattens 
out at its lowest value. Both the minimum and maximum 
temperatures follow the same behavior. The distortion and 
average temperature rise displayed in Fig. 16 show the same 
patterns. Note that for this configuration of model and alti- 
tude, the jets will converge at the ground for a thrust splay 
angle of 22 deg. 

Splaying Forward Jets Only 

In the second variation, only the forward jets varied in 
thrust splay angle. The rear lift jets maintained a 0-deg splay 
angle. Again, the thrust splay angle varied from 0 to 45 deg 
(5 = 0-45 deg) in 5-deg increments at the same altitude (H 
= 4 D s ) and head wind speed (LL = 0.03 V,), just as in the 
all jets splayed variation. The calculated results are extremely 
similar to the results from splaying all the jets. As a matter 
of fact, a plot of the temperature contours in the same two 
planes as in Fig. 14 would appear almost identical to Fig. 14 




Fig. 16 Distortion and average rise vs splay angle on all jets. 
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Fig. 17 plot of engine face temperatures vs splay angle on forward 
jets. 
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Fig. 18 Distortion and average rise vs splay angle on forward jets. 


67 


forward of the jets. (Such a plot does exist in the paper, but 
it is not here for brevity.) 

The ingestion effects of thrust splay angle for the forward 
jets can be seen in Fig. 17. The average temperature first rises 
and then drops until a splay angle of about 20 deg, where it 
rises very slightly. A local maximum exists at 25-deg splay 
angle before the average temperature drops and flattens out 
at its lowest value. Both the minimum and maximum tem- 
peratures follow a similar behavior. The distortion and av- 
erage temperature rise displayed in Fig. 18 show the same 
patterns, just as when all the jets are splayed. Overall, splay- 
ing the forward jets alone gives the same effects on hot gas 
ingestion as splaying all the jets. 


Conclusions and Remarks 

In summary, the average of the engine face temperatures 
decreases with increasing height, but is relatively unaffected 
by head wind speed. Engine face temperature distortion also 
decreases with increasing height, but increases with head wind 
speed until the forward vortex is behind the inlet face. The 
head wind speed variation reveals that (for a constant height) 
the engine face temperature is dominated by near-field inges- 
tion effects. 

As for the thrust splay angle variation, splaying the jets 
inward (for a constant height and head wind speed) first causes 
a rise and then a rapid decrease in hot gas ingestion, although 
a local maximum exists when the jets converge at the ground 
plane. Also, splaying the forward jets alone, instead of all 
the jets, gives almost the same reduction in hot gas ingestion 
without as large a thrust penalty. 

A comparison of the H = 4D ; , U — 0.03Vj cases in the 
height parameter variation and the head wind speed param- 
eter variation show the effects of the calculation domain on 
the flowfield. The vertical walls of the narrower calculation 
domain definitely affect the ground vortices, but no differ- 
ences exist in the temperatures reaching the engine. If the 
overall flowfield is of primary interest, then the tunnel walls 
would have to be farther from the aircraft model. 

This study did not address the importance of the aircraft 
geometry (fuselage, wings, tails, etc.) in relation to the flow- 
field. Only one aircraft model was used, and it was quite 
simplistic. . 

The last conclusion from this study concerns the practicality 
of using an efficient CFD code for parameter variation studies. 
The turnaround time on a Cray-2 supercomputer and state- 
of-the-art workstations allows quick parameter changes. Typ- 
ically, a Cray-2 supercomputer solved the flowfield in about 
an hour with a turnaround time of a day. A dedicated IBM 
RS-6000 workstation can solve the flowfield in about 6 h and 
can actually give shorter turnaround than the shared super- 
computer. 
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Abstract 

The development, validation, and application of a general 
purpose multigrid solution algorithm and computer program for 
the computation of elliptic flows in complex geometries is 
presented. This computer program combines several desirable 
features including a curvilinear coordinate system, collocated 
arrangement of the variables, and Full Multi-Grid / Full 
Approximation Scheme (FMG/FAS). Provisions are made for 
the inclusion of embedded obstacles and baffles inside the flow 
domain. The momentum and continuity equations are solved in 
a decoupled manner and a pressure correction equation is used 
to update the pressures such that the fluxes at the cell faces 
satisfy local mass continuity. Despite the computational 
overhead required in the restriction and prolongation phases of 
the multigrid cycling, the superior convergence results in 
reduced overall CPU time. The numerical scheme and selected 
results of several validation flows are presented. Finally, the 
procedure is applied to study the flowfields in a side-inlet dump 
combustor and twin jet impingement from a simulated aircraft 
fuselage. 


Introduction 

In recent years, multigrid methods have been shown to 
effectively improve the speed of many internal and external 
fluid flow calculations 1 * 4 . * Multigrid methods have been 
applied to both incompressible and compressible flows and 
unstructured grids 5 . Despite these advances, there have not 
been many applications of multi grid methods to complex 
practical internal flows that include features such as embedded 
obstacles, mass injections, combustion and complex 
geometries. 

In the present paper, we describe our recent work towards the 
development of a finite volume multi grid algorithm as well as a 
computer program that can address several practical flows 
encountered in the propulsion industry. The solution algorithm 
is based on a collocated arrangement of pressures and velocities 
and uses multigrid cycling to efficiently converge the low 
frequency errors. A combined flrst/second order difference 
scheme is used to discretize the equations. The steady-state 
mean flow equations along with a two-equation k-e model of 
turbulence are solved for analyzing turbulent flows. Further, 
provisions are made to include any arbitrary number of 
obstacles and baffles to simulate internal flow obstructions. 

The concept of multigrid methods is to use a sequence of coarse 
and fine finite-difference grids on which the solution is cycled 
continuously to remove the low frequency content of the 
solution error. The current multigrid cycling uses a V-cycle for 
transferring the residuals and prolongating the corrections from 
the coarse grids. In the collocated arrangement, the velocities 
and pressure are located at the cell centers, but the fluxes are 
evaluated at the cell faces. For incompressible flows, this 
evaluation of fluxes must be done carefully in order to avoid a 
checker-board split of the pressure field. The present scheme 
uses a momentum interpolation scheme similar to the practices 
of Rhie and Chow 6 , Peric et al. 7 , and Barcus et al. 4 
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In the following sections, a description of the governing 
equations, the numerical procedure, the systematic validation of 
the code, and its application to two problems of interest to the 
aero-propulsion community are presented. 

Governing Fluid Flow Equations 

All the relevant partial-differential equations can be represented 
by a single equation for the transport of a general scalar variable 
<J>. The governing equation for <f> in an arbitrary coordinate 
system (q,q,Q can be written in strong conservative form as 

^ (pu<t>) + (pv« + |r (pw<0 = 

dq dr| dq 

^[j(qn + <U2 <t>Ti + <U3 < t>£)] + 

^[j (q21 H + <122 <t>Ti + q23 ♦(;)] + 


^[j (q31 H + q 32 *T1 + q33 ♦{;)] 


Sft.il. CW (1) 

where U, V and W are the contravariant components of the 
velocity vector. These components are related to the Cartesian 
velocities by the relations 

[ li T 3 = J [ A ] [ jiT ] (2) 
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where 

Jl T = (u, V, w) 

The metrics in these expressions are given by 


(3) 
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, 2 2 2 ,.2 

fl22 = ( T 1x + ^y + T1 z > J (8) 

etc. 

For the turbulence transport equations, additional production 
and dissipation source terms in the equations for k and e are to 
be included. These terms are transformed to the curvilinear 
coordinate system by the chain rule of differentiation. 

Numerical Procedure 

For a collocated scheme, the mass fluxes at the cell faces and 
the velocities at the cell centers are separately calculated and 
stored. The algorithm currently employs the popular hybrid 
differencing for evaluating the interface values. The discrete 
equation can be written in the conventional form 

ap<t>p = £ a n fc> <(>nb + ^ W 

nb 

where <{>nb are the neighbor values of <t>p on a seven point 
stencil. SV is the corresponding total source term. 

The key feature in collocated schemes for incompressible flows 
is the momentum interpolation of two cell-centered Cartesian 
velocities. The central idea of momentum interpolation^ is to 
alter the expression for the net dp/9^ at the (i+1/2) cell face with 
a staggered pressure difference. For the V and W fluxes, the 
effective dp/orj and dp/dt, are similarly modified by staggered 
pressure differences. Such a momentum interpolation then 
provides a well-connected pressure field. 

Consistent evaluation of the volume fluxes (U, V, W) on a 
coarse grid and prolongation of the associated corrections to the 
finer grids are observed to be the most important issues in the 
present collocated multigrid procedure. An inconsistent 
procedure resulting in limit cycling of the mass residuals is 
observed if the coarse grid fluxes at the cell faces are not 
properly updated after the momentum equations are solved. Let 
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be the restriction of the fluxes and the Cartesian velocities from 
a fine grid h to a coarse grid H, respectively, during the 
multigrid cycle. The solution of the coarse grid momentum 
equations gives a velocity field that satisfies 

L H u H = f h + R h (12) 

where R^ is the residual of the coarse grid equation given by 

gH = I« Rh (13) 

Let 

5uH= u H-SH (14) 

A consistent update of the coarse grid flux is obtained by 
evaluating 12^ as 

jjH - £)H + J • [ A ] [ 5 ijH ] T (15) 


where [ A ] is the matrix of coordinate transformation on the 
coarse grid. This coarse grid flux, however, does not satisfy 
the local mass continuity equation, thus requiring the restricted 
pressure field 

pH = l”ph (16) 

to be corrected by a pressure-correction equation. The coarse 
grid pressure-correction equation is derived in an identical 
manner to the fine grid equation by implying a staggered 
location of the cel! face fluxes. 

It is necessary to point out that the direct evaluation of the 
coarse grid fluxes from the newly computed Cartesian 
velocities leads to an inconsistent formulation with the residual 
mass source reaching a fixed non-zero value. Thus, the 
seemingly straight-forward update procedure 

U H = J - [ A ] [ U H 1 T 07) 

on a coarse grid during restriction phase does not lead to a 
consistent formulation. The derivation of the pressure 
correction equation follows the perturbation concept introduced 
by Patankar and Spalding 8 and used widely. The complete 
details of the solution algorithm are given in Smith and Vanka 9 . 

Results 

The solution algorithm and its multigrid efficiency have been 
tested in a number of model problems. A detailed description 
of the systematic validation is precluded by space limitations; 
therefore, only three selected calculations are presented. 

Laminar Lid-Drive n Cube 

To demonstrate the efficacy of using multigrid methods on 
elliptic problems, the laminar flow in a lid-driven cube was 
simulated for a Reynolds number of 100 based on lid velocity 
and side dimension. A 32 3 cell solution was selected. This 
problem was solved using identical initial guesses, relaxation 
factors, convergence criteria, and other numerical parameters in 
both a single-grid mode and in a three-grid mode. Figure 1 
presents the convergence history for both a single grid solution 



Fig. 1 Convergence history for laminar lid-driven cube 
calculation. 


and a multi grid solution using a standard (1,1,1) three-level V- 
cycle. The multigrid solution converges at a log-linear rate on 
all three grid levels, while, in comparison, the single grid 
solution has an oscillatory and sluggish rate of convergence. 
The multigrid solution required only 20 equivalent iterations to 
converge to a mass residual of 10~->; however, the single-grid 
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solutionrcquired 1 22 iterations to reach the same convergence 
level The multigrid solution represents a speedup of 3.75 
(based on CPU timings) over the single-grid solution. Further 
speedup can be expected when finer grids are considered. 

Laminar Pjpp fifnfl 

As a stringent test to the coordinate transformations, the laminar 
flow in a 90 degree strongly curved pipe was analyzed. The 
now in strongly curved ducts is elliptic because of the strong 
radial pressure variations and possible axial flow 
recirculation 10 . The circular pipe was transformed to a square 
through an elliptic grid generation procedure 1 The curvature 
of the duct is also represented by a grid transformation to a 
straight duct. This flow has been previously computed by 
other researchers but on coarser grids and using single grid 
techniques. 

The calculation was performed for the geometry described by 
Enayet et al.12 In these experiments, the radius ratio of the 
bend was 5.6 and the Reynolds number was 500 
corresponding to a Dean number of 211. The computational 
domain consisted of a short upstream tangent, the bend itself, 
and a downstream tangent. A representative computational grid 
is shown in figure 2. The advantage of the present coordinate 
system versus the conventional (r, 0, $) system used in earlier 
works is that the current system distributes the mesh points 
more uniformly in the cross-section and avoids numerical 
difficulties caused by large cell aspect ratios at the center. The 
finest grid used in the present study consisted of 120 x 64 x 32 
cells in the streamwise and cross-sectional directions, 
respectively. In order to accurately prescribe the thick inlet 
boundary layer, a fourth-order polynomial was fitted to the 
experimental streamwise velocity data at 0.58 diameters 
upstream of the bend inlet plane and prescribed as the inlet 
velocity to the computations. The remaining boundary 
conditions were no-slip at the walls, symmetry conditions at the 
center plane, and a zero-derivative condition at the exit. 

shows the rate of convergence for a three-level 120 x 
o4 x 32 computation. The coarsest grid converges slowly due 
to the poor initial guess and is not shown. However, the two 
finer grids converge to a good accuracy at essentially the same 

v/? 4 D A ^ roximale ^ ^ m * nutcs wc re required on a CRAY- 
Y/MP computer to complete the computation. Sequential grid 
refinement demonstrated that 120 x 64 x 32 cells provide an 
adequate gnd-independent solution. Figure 4 presents a 




Fig. 3 Convergence history for laminar curved pipe 
calculation. 



Fig. 4 Calculated and experimental streamwise velocity 
profiles at the 60° plane in the laminar curved pipe 
( radial coordinate versus non-dimensional 
streamwise velocity). 

comparison of calculated values over a range of grid 
distributions and the experimental data of Enayet et al. ^ at the 
W degree plane. The agreement is very good with the high 
peaks being extremely well predicted. Agreement at other bend 
angles is observed to be also equally good. Figure 5 shows 
cross-stream flow patterns at 60 degrees indicating the 
production of streamwise vorticity by the bend. 



Fig. 2 Representative computational grid for curved circular 
pipe calculation. 
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Fig. 5 Calculated secondary flow at the 60° plane in laminar 
curved pipe calculation. 



Turbulent Flow in a C urved Square Duct 

As a final validation case, calculations of the turbulent flow in a 
strongly curved square duct were performed. The 
computations were performed for the geometry described by 
Taylor et al. 13 including a 2.3 radius ratio at a Reynolds 
number of 40,000. Again, upstream and downstream tangent 
sections were included in the computational domain. The 
turbulence modeling is based on the standard k-E model 
without modifications to account for the curvature-induced 
anisotropy effects. Admittedly, this will result in some 
discrepancies between measurements and calculations. The 
solution of the k-e equations in conjunction with the momentum 
equations is carried out in a sequential manner, but the k-E 
equations are solved only on the locally finest grid in the Full 
Multigrid cycle and the turbulence viscosity field is restricted 
from the fine to coarse grids. Wall functions 14 arc imposed on 
the finest grid and the implied wall viscosity is restricted for use 
in the coarse grid momentum equations. The finest grid in the 
present study utilized 120 x 64 x 32 cells in the streamwise, 
spanwisc, and transverse directions, respectively. 

Figure 6 shows the convergence of the algorithm on the two 
finest grids. It is seen that the convergence is not as fast as for 
the laminar case; however, the number of iterations required is 
still small and independent of the grid density. The solution 
required approximately 20 minutes on a CRAY-Y/MP 
computer. The slower convergence for turbulent flows in 
comparison with laminar flows is due to the coupling between 
the momentum and turbulence equations. It may be possible to 
accelerate this convergence slightly by better optimization of the 
relaxation factors. 



Fig. 6 Convergence history for the turbulent curved square 
duct calculation. 


Figure 7 shows comparisons between calculated and 
experimental streamwise velocity profiles at angular positions 
of 30 degrees and 77.5 degrees in the bend. The calculations 
show satisfactory agreement with measured values and are also 
in good agreement with the calculations of Kreskovsky et al. 
who used a similar model for the turbulent closure. However, 
the agreement with experimental data can be improved with 
more advanced. turbulence models such as those based on the 
solution of the Reynolds stress equations. 

Side-Inlet Ducted Rocket 

As a demonstration of the computer program, the flowfield in a 
side-inlet ducted rocket was studied. The investigation was 
limited to an isothermal water tunnel simulation as reported by 
Stull et al. 16 This turbulent flowfield possesses several 
complex features including a streamwise recirculation region 
within the dome, impinging jets within the inlet region, and a 
pair of cross-stream vortices downstream of the inlet. 




d' b 0 

(a) (b) 


Pig. 7 Calculated and experimental streamwise velocity 
profiles along a mid-span line in turbulent curved 
square duct at a) 30 degrees and b) 77.5 degrees 
( transverse coordinate versus non-dimensional 
streamwise velocity). 


The computations were carried out over a grid of 96 x 48 x 24 
cells in the streamwise, transverse, and spanwise directions, 
respectively. Figure 8 illustrates a representative computational 
grid in the cross-sectional plane. Boundary conditions included 
no-slip walls for the combustor lining and dome plate, 
symmetry for the mid-span plane, inflow segment for the inlet 
duct, and zero-derivative outflow for the exit plane. The sice- 
inlet arms were modelled as an inlet boundary condition with a 
plug velocity of 2.67 m/s corresponding to a inlet duct 
Reynolds number of 2.0x1 0 5 . Three grid levels were used and 
the solution was iterated until a normalized mass residual of 
1 OxlO- 3 was obtained. Figure 9 illustrates the good 
convergence rate on all grids. The slow convergence towards 
the end on the third grid is believed to be a result of the 
coupling between the turbulence and flow equations as 
described previously. Approximately 20 minutes on a CRAY- 
Y/MP computer were required to complete the computation. 



Fig. 8 Representative cross-sectional computational grid for 
ducted rocket calculation. 
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Iteration 

Fig. 9 Convergence history for ducted rocket calculation. 

Figure 10 presents the computed velocity vectors in the 
symmetry plane of the ducted rocket. The entrainment by the 
inlet jets is responsible for driving the dome recirculation 
pattern shown. Figure 11 shows computed cross-stream flow 
patterns at two axial locations. Downstream of the inlet region, 
the flow is seen to possess a cross-stream vonex which results 
in helical pathlines. 

The authors are not aware of any existing quantitative 
experimental data to which the present calculations can be 
directly compared; therefore, only qualitative observations are 
reported. The computed flow patterns agree well with flow 
visualization tests 16 except in the dome region where the flow 
was observed to be unsteady. The present calculations solve 
the steady-flow equations and, therefore, are not expected to 
give good agreement in this region. The flowfield downstream 
of the inlet is well characterized by the computations. 

STOVL Hot Gas Ingestion Study 

The eventual application of this computer program is intended 
to be the study of ingestion of hot gases by the inlets of Short 
Take-off and Vertical Landing (STOVL) aircraft in ground 
proximity. A number of interesting fluid dynamic effects have 
been identified when the lift jets of STOVL aircraft impinge on 
the ground surface 17 * 18 . First, as the lift jets expand, they 
entrain the surrounding fluid causing a negative pressure 
underneath the fuselage and a loss in lift. As the jets impinge 
on the ground and spread radially outwards, the wall jets 




lug. 1 1 Computed cross-secdonal velocity vectors in the 

ducted rocket a) immediately downstream of inlet and 
b) at 1.0 diameter downstream of inlet. 


Fig. 10 Computed velocity vectors in the symmetry plane of 
the ducted rocket. 
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further entrain external fluid and increase the loss in lift. In a 
multiple-jet configuration, these wall jets-collide with each other 
and turn upwards to form an upwash fountain. This fountain 
flow has two main effects on the STOVL aircraft dynamics. 
First, an increase in lift force is caused when the fountain 
impinges on the aircraft fuselage. The'rccovery in lift is a 
positive effect of the upwash flow. However, this impinging 
fluid can also flow along the fuselage surface and eventually 
make its way into the engine inlets. Because the temperature of 
the fountain flow is much hotter than the ambient air, its 
ingestion by the engine can reduce the power and cause thermal 
stresses in the components. In addition to the fountain flow, 
another mechanism for the hot gas ingesoon results from the 
interaction of the forward moving wall jet with the headwind. 
When the headwind and the wall jet collide, a stagnation region 
is formed and the wall jet is turned into a ground vortex. This 
ground vortex can subsequently be ingested by the engine 
resulting in a further loss in power. 

Several numerical studies of twin- and single-jet impingement 
have been previously reported 18 * 19 . Recently VanOverbeke 
and Holdeman 20 studied the hot gas ingestion process in a 
simulated geometrical configuration of the STOVL aircraft. 
The study was restricted to Cartesian geometries but 
demonstrated the potential of computanonal fluid dynamics to 
this practical flow problem. This study was followed by Tafti 
and Vanka 21 who applied a multigrid solution procedure and 
demonstrated significant reductions in computer times over a 
single grid procedure such as that used by VanOverbeke and 
Holdeman. In this section, we present results of a preliminary 
calculation of the impingement of twin jets from a simulated 
fuselage represented by a curvilinear grid. 

The simulated aircraft fuselage was a truncated NACA 0012 
airfoil at 10 degrees angle-of-attack embedded within the 
computational domain. Figure 12 shows the grid in the 
impingement region. The calculation domain extended far 
upstream but is not shown in this figure. The lift jets were 
specified as 300 m/s at 1000 K and the headwind was 
prescribed as 3% of the jet velocity at an ambient temperature of 



300 K. The flowfield was considered to be turbulent. A 104 x 
32 x 32 grid distribution with two grid levels were used in this 
calculation requiring approximately 150 iterations to converge 
to a normalized mass residual of one percent. The above 
calculation required approximately one hour on the CRAY-2 at 
NAS of NASA-Ames Facility. 

Figure 13 shows a selected velocity field in a plane through the 
jet centers. The impingement of the jets on the ground and the 
formation of the wall jets is evident in this plot. The headwind 
is seen to create a stagnant zone ahead of the fuselage and the 
jet flow is turned back by the headwind. The interaction of the 
hot jets with the cooler ambient is illustrated in Figure 14 in two 
different planes. In figure 14-a, the temperature distribution in 
a plane through the jet centers is shown- In the present 
calculation, engine suction was not considered; therefore, the 
temperatures at the engine face could not be estimated. This 
feature will be exercised in future calculations. Figure 14-b 
shows the temperature footprint of the lift jets on the ground 
plane. The interaction of the wall jets with the headwind can 
also be seen in this figure. 


Summary 


A multigrid based computational algorithm and code have been 
developed and validated for elliptic complex flows. The 
algorithm solves the momentum equations for the Cartesian 
velocities as the dependent variables and stores all the variables 
at the cell centers. A consistent multigrid formulation in 
curvilinear geometries has been developed and tested. The 
rapid convergence of the algorithm has been demonstrated in a 
variety of complex flows. 
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Fig. 1 4 Contours of temperature in a) the jet central plane and 
b) ground plane for twin-jet impingement flow. 
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Abstract 

A multigrid method for internal flows in complex geometries 
based on multiple velocity grid staggering has been 
described. The numerical method has been tested for 
laminar flow in a curved pipe and for laminar and turbulent 
flow in the dilution zone of a gas turbine combustion 
chamber. The convergence characteristics for the method 
have been compared to the convergence obtained for use of 
traditional grid staggering. It has been observed that the use 
of the multiple velocity grid staggering arrangement is 
superior to traditional grid staggering only when a ninety 
degree bend exists in the problem geometry. For the 
simulation of flow in the dilution zone of a gas turbine 
combustion chamber, it was observed that both gnd 
staggering arrangements required approximately the same 
number of iterations for convergence. In this case, the 
traditional grid staggering is to be favored due to its lower 
work count per iteration. 

I. Introduction 

Multigrid methods have long been advocated for 
accelerating the convergence of iterative schemes used for 
the numerical solution of elliptic partial -differential equations 
(Brandt, 1977; Brandt, 1980). Multigrid methods rely on 
the simultaneous use of several grid levels on which the 
governing equations are discretized and iteratively solved. 
The principal argument for multiple grid levels is that the 
various frequencies in the error spectrum are efficiently 
resolved at the corresponding high frequency rate of the 
iterative scheme. Thus, a stria O(n) rate of work increase 
can be achieved. The theory of multigrid methods is well 
established for model linear elliptic partial-differential 
equations (Brandt, 1980). Also, multigrid methods have 
been applied to a large number of linear and nonlinear 
partial-differential equations including the Navier- Stokes 

equations. . 

The application of the multigrid technique to accelerate 
the calculation of practical engineering fluid flows is not well 
established. Although a number of earlier works (Vanka, 
1986a; Rayner, 1991; Rubini et al., 1992; Shyy ct al., 1993; 
etc.) have solved the Navier-Stokes equations for mode 
flow geometries, the application of multigrids to practical 
flows that necessitate the solution of a much larger set of 
coupled equations in complex geometric has not been fully 
investigated. Cunendy, computational fluid dynamics plays 
an important role in industry in the evaluation of proposed as 
well as existing designs of engineering components, today, 
commercial software packages are in the position of 
providing answers to several questions arising outo ^ n ^ 
designs. Therefore development of efficient computational 
algorithms utilising the multigrid concept for pracucai 
engineering flows is of much significance. 

Ln several of our previous works (Vanka, 1986b, Joslu 
and Vanka, 1991), we had developed mulngnd based 
algorithms for flows in complex geometries with inclusion 
of models for turbulence and combustion. While these 
works were satisfaaory in their own respect, the issue of the 


* Teaching Assistant, Member AIAA 

** Research Assistant 

*** Professor, Dept, of Mech. and Ind. Eng. 
Senior Member, AIAA 


80 


treatment of complex geometries still remains to be properly 
resolved. Our interest is primarily with methods involving 
structured grids such as those generawl from the solution of 
elliptic partial differential equations (Thompson et al., lysoj, 
rather than unstructured grids such as those used in the 
finite-element methodology. The key issues in our previous 
works have always been the selection of the proper set o 
dependent variables for the momentum equations and the 
layout of the velocities with respect to the lwanonsof the 
pressures (Joshi and Vanka, 1991; Smith et aL, 1993). This 
issue is of much more significance to incompressible flows 
as it is possible to generate a checker-board split m the 
pressure field, if the pressure gradient terms in the 
momentum equations are not properly discretized. The 
various schemes have ranged from those using the Cartesian 
velocities located on the cell faces and pressures at the cell 
centers (Vanka et al., 1989) to the use of a collocatea 
arrangement of velocities and pressures at the cell centers 
(Smith et al., 1993). While these practices have had 
ad rq yare success in solving complex flows, they have also 
resulted in certain difficulties. 

The location of single Cartesian velocity components on 
appropriate cell faces has been pursued by Vanka et al. 
(1989) and Shyy and Vu (1991). Its mulngnd performance 
has been discussed in Vanka (1987) in the context of a 
block-implicit solution algorithm. This algorithm however 
has been observed to have convergence difficulties when the 
angle between the grid lines is beyond a certain value (~ 40 
degrees). In cases when the flow geometry turns ninety 
degrees or more, the layout can result in a staggered 
arrangement that is contrary to the desired objective. The 
use of collocated Cartesian velocities at the cell centers 
resolves the difficulty arising from arbitrary gnd skewness 
and can be incorporated into a multigrid procedure with 
sequential solution of the momentum and continuity 
equations (Miller and Schmidt, 1988; Hottmann et al., 1990, 
Smith et al., 1993). However, in selected problems, we 
have observed that this procedure can result in 
inconsistencies during the multigrid restriction and 
prolongation operations. This causes the overall algorithm to 
converge to a limit residual (-1%). Further, the restriction 
and prolongation operators must always be consistent with 
the momentum interpolation that is used to avoid the 
checker-board decoupling of the pressure. The use of ceU 
face mass fluxes or grid aligned velocities (Karkj and 
Patankar, 1988; Joshi and Vanka, 1991) as dependent 
variables can be another alternative, but both practices lead to 
complex forms of curvature terms which again are difficult 
to treat in a multigrid context. 

In the present work, we have considered an extension of 
a practice originally proposed by Maliska and Raitiiby 
(1984). The proposed scheme locates all components ot the 
Cartesian velocity vector on each cell face and solves the 

respective momentum equations. Thus in a two-dimensional 

problem, two Cartesian velocities are located on each cell 
face. Four momentum equations are solved tor tne 
components of velocity on the £ and 7 } cell faces and are 
directly used in the continuity equation. This arrangement 
avoids the difficulties associated with grid turning and does 
not require momentum interpolation as used in the collocatea 
scheme. However, its disadvantage is the additional work 
and storage for the multiple momentum equations. This 
work at first appears to be twice (two dimensional problems) 
or three times (three dimensional problems) the work of a 
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traditional single component algorithm. However, if we 
consider the fact that a common set of finite-difference 
equations can be assembled for all equations on a given cell 
face and that the momentum equations require fewer sweeps 
(in the context of the SIMPLE algorithm (Patankar and 

Spalding, 1972)) than the pressure and k-e equations, it can 
be shown that the overall increase in work is at most a factor 
of two in three dimensional problems and less in two 
dimensional problems. Similarly, the storage increase is 
also not large in a relative sense. 

In the present paper, we describe our work relating to the 
development and assessment of a multigrid algorithm for 
complex internal flows that solves multiple Cartesian 
velocity components on the cell faces. The algorithm uses a 
segregated iterative scheme for resolving the pressure- 
velocity coupling. The multigrid procedure uses a FAS- 
FMG strategy and a V-cycle for restricting and prolongating 
between grids. The discretization uses either the traditional 
hybrid differencing or a third-order flux limited scheme 
(Leonard and Mokhtari, 1990). The algorithm has been 
incorporated such that with one index it is possible to 
retrieve the single component algorithm appropriate for 
simple geometries. We have applied the procedure to two 
model flows: a) laminar flow in a curved pipe and b) the 
flow in a model annular combustor. We have solved the two 
problems in both single grid and multigrid modes and with 
single and multiple Cartesian velocity procedures. As 
expected, the multigrid algorithm converges faster than the 
single grid algorithm, although the benefits are not as great 
as would be desired due to the relatively coarse 
demonstration grids considered. A two equation turbulence 

model (k-e) has also been included in the solution procedure 
but currently the performance of the overall algorithm for 
turbulent flows is in the preliminary stage. 

In what follows, the discretization procedures and the 
governing equations are described in section 2. The solution 
procedure is treated in section 3. Sections 4 and 5 present 
the results for the two problems considered in this study. 

IL Governing Equations and Discretization 

Governing Equations and Boundary Conditions 

The conservation laws can be written in terms of a 
general equation containing convection, diffusion, and 
source terms. For a general variable <p we solve the 
following equation: 


-^(pu<t>) + jj(P v *) + = 
d r r* 
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where the metrics of the coordinate transformation and 
contravariant velocities are given by the following 
expressions: 
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In the present investigation, turbulence closure has been 

accomplished by using the k-e model proposed by Launder 
and Spalding (1974). 

The numerical method can accommodate a variety of 
boundary conditions including Dirichlet, Neumann, and 
outflow boundary conditions. For the problems addressed 
in this study, the values of velocity, temperature, density, 
and turbulence variables are specified at inlet boundaries. A 
Neumann boundary condition at the outflow plane requires 
that the normal derivative of all velocity components, 

temperature, k, and £ be equal to zero. Wall boundaries are 
treated as no slip, no penetration surfaces. The wall 
boundaries arc treated as either constant temperature or 
adiabatic surfaces for the solution of the energy equation. 
For the simulation of turbulent flow, the wall functions 
proposed by Launder and Spalding (1974) have been used 
for the three velocity components and temperature. The wall 
functions have been implemented by writing the effective 
viscosity in terms of the wall shear stress, where the wall 
shear stress is obtained from the turbulence equilibrium 
assumption. 

Discretization 

A multiple velocity grid staggering scheme has been used 
in discretizing the governing equations. All three Cartesian 
velocity components are stored on the cell faces of a given 
control volume (see Figure 1). The pressure, density, and 
other scalar variables are located at the cell centers. The grid 
staggering is an extension of the scheme proposed by 
Maliska and Raithby (1984). For a given cell, nine 
Cartesian velocity components are solved. The multiple 
velocity grid staggering has been utilized because it does not 
exhibit pressure- velocity decoupling when ninety degree 
bends occur in the geometry. Furthermore, the velocities 
which arc solved in the momentum equations are the same 
velocities updated in the pressure correction equation, unlike 
the situation for collocated grids. Thus, the same velocities 
will satisfy both the momentum and continuity equations. 

The governing equations are discretized by a finite 
volume method. The integration of the partial differential 
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equation over a control volume results in the following 
equation: 


limited through use of normalized variables, where a 
normalized variable is defined by the following: 


[ipu*)r ~ (pW) f - ] + [(P^)„- - (P v *\- } 


[(n/-V„^) r -{rr'q u ^) v 
[(r/-‘^) 5 . -(rr^ A l 
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where S^ D refers to the mixed derivative terms in curvilinear 
coordinates. This equadon is discretized by specifying the 
relationship between convective and diffusive cell face fluxes 
and the surrounding volume averaged quantities. The above 
equation is first linearized by treating the cell face mass 
fluxes and exchange coefficients explicitly. These cell face 
quantities are approximated by using linear interpolation. 
The cell face values of the dependent variables and diffusive 
fluxes are then discretized and treated implicitly. The 
diffusive flux at a control volume face is discretized by a 
second order central difference approximation: 


d0 = <P E -<Pp 

aft *5 


The dependent variable at a cell face is described by the 
hybrid differencing procedure (Spalding, 1972). As the 
discretized equations using this differencing are well known 
and thoroughly documented they will not be repeated. The 
results of calculations using a higher order discretization are 
also given and the discretized equations using the higher 
order flux function will be presented. 



In the higher order differencing scheme, the cell face 
value of a given dependent variable is related to the 
surrounding volume averaged quantities by a flux limited 
form of quadratic interpolation for convective kinematics, 
QUICK (Leonard, 1979). In order to rigorously treat a cell 
face quantity using this method, a multi-dimensional 
interpolation operator should be used (Leonard, 1988). In 
the present study, the one dimensional form of quadratic 
interpolation has been implemented. For the cell face value 

of a variable 0 y the following expression is used: 

< t>/ = ^ < t>D + ^ < t>c-J l Pu (5) 

Where 0 Dy 0 uy and 0 C are the downstream, upstream, and 
central values of the dependent variable. The cell face flux is 



At a given cell face, the downstream, upstream, and central 
values of the dependent variable are identified. The 

normalized variable associated with the center quantity, 0 C , 
is used to assess the local behaviour of the solution of 0 . 
Depending on the local variation, the value of the cell face 

normalized variable, 0 fy is determined. In this study, the 
universal limiter of Leonard and Mokhtari (Leonard and 
Mokhtari, 1990) has been utilized. The value taken by <p f is 
determined from 0 C in the following manner 
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The functional dependance of 0 f on 0 C is depicted in Figure 
2. In order to implement the flux limited interpolation 
implicitly, the downwind weighting factor is defined: 

DWF = ~ <t> - = (12) 

Qd-Qc l ~<t>c 

Having determined 0 f at a given cell face from equations 7- 
1 1, the downwind weighting factor can be calculated and the 
cell face value of the dependent variable written as follows: 

0 f = DWF 0 D +(l- DWF ) - 0 C (13) 



Figure 2: The Universal limiter constraints imposed on a cell 
face variable as shown on the Normalized Variable Diagram 
(NVD). 

A deferred correction procedure has been used to 
implement the flux limited QUICK formulation. The 
convective cell face flux is written in terms of the first order 
upwind value plus the difference between the higher order 
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and first order approximations. For a q -plus cell face flux, 
the following expression is obtained: 


(pU0) r = [<t>p max[(p£/) r ,o]- <p E ■ max[-(pt/) r ,o]| 

| DWF r - <t> e + (l - DWF r )<pp\ ■ max[(pt/) r ,o] - 

\DWF V <t> P + ( 1 - DWF t . )* e ] • max[-(pt/) r ,oJ 

-\<p p ■ max[(p(/),. ,0| - <!>e ' max[-(pt/)^. ,ojj 
1 L (14) 


The first order upwind term is treated implicitly and the 
difference term is treated explicitly as a source. Such a 
procedure maintains diagonal dominance in the system of 
linear algebraic equations and ensures thatthe coefficients 
are always positive (Hayase et al., 1992). The final form o 
the discretized equation for the variable p is thus given by 

the following: , , . 
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alternating line Gauss-Seidel procedure were performed for 
each velocity and scalar variable while ten iterations wer 
performed on the pressure correction equation. The six 
steps given above constitute one "outer iteration, l-. 
iterative update of all flow variables. A numerical me 
employing only a single grid would simply continue these 
outer iterations until the residuals in the governing equations 
were reduced below a prescribed level. In the present 
multigrid method, the outer iteranons are performed on a 
given grid prior to restriction or prolongation to “°^gnd. 
It is important to note that during the execution of a V-cycle, 
the scalar equations are only solved on the locally finest gnd. 
The values of T, k, and e on the finest grid are held fixed on 
the coarse grids during the V-cycle. The details pertaining to 
the pressure correction equation are described, below. 

Pressure Correction Equation 

The pressure correction equation used in the present 
study is obtained by writing expressions for the mass fluxes 
associated with a given cell in terms of an ininal esnrnate 
plus a correction. The integral form of the conservation of 
mass for a control volume is given by the following. 


[(pt/) r -(p(/) r ] + 

[(pv),.-(pv) r ]+ 

[(pW) { .-(pW) 5 -] = ° 


(16) 


The relation between flux corrections and pressure 
corrections is obtained from the momentum equanons. 
Consider the mass flux at the £-plus cell face. The 
correction to the mass flux is written in terms of velocity 
corrections 

pt/ r =p(y„z { - V-) r -(4 + “r) 

+p{ x f z v - -V«) r ' ( v r + v r ) (17) 
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III. Solution Procedure 

The coupled set of nonlinear algebraic equations which 
results from discretization is solved using a multignd 
solution procedure. The relaxation on any given gnd is 
accomplished by using a segregated solution technique 
(Patankar and Spalding, 1972). A FMG-FAS multignd 
procedure is used to accelerate the rate of convergence of the 
basic relaxation operator. 


Relaxation Operator 

In the segregated relaxation method, each equation is 
solved in a sequential fashion. On a given gnd, the 
following procedure is used: 

1) Solve for all Cartesian velocities at £ faces. 

2) Solve for all Cartesian velocities at tj faces. 

3) Solve for all Cartesian velocities at q faces. 

4) Solve for the pressure correction, P'. 

5) Update Cartesian velocities and P based 

on the solution for P\ 

6) If on the locally finest grid, solve for 

scalar variables (T, k, e) 

The iterative process for any particular flow variable is not 
carried to convergence since subsequent iterations on the 
other flow variables will change the coefficients of the 
discrete equation. In general, three to five iteranons of an 


Since all three Cartesian velocities are stored at any cell face, 
the relations for u\ v\ and w’ at the 4 -plus cell face are 
obtained from the momentum equations for the velocites u, 
v, and w at that cell face. The following expression for U v 

is thus obtained: 

tf { .=(v 5 -y< z O r V 

— w f 
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should be noted that for simplicity, the cross derivative 

-ms in pressure {dP'/dn, 3P7*) arc neglccted , 
jwever, when the included angle between *jy <wo grid 
ies is less than thirty degrees, the use of the cross 
rivative terms may be needed to accelmte die 
eric 1990). The expressions for V and W are 
SI k similar logic. The final pressure correction equation 
suiting from the perturbed continuity equation can be 
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written as follows: 

a p P’ p = a E P' E +a v P' w + a N P^ + a s P' s + a H P' H + a L P' L - m„ 
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approximarions on grid (k-1) to the correct values on grid 
(k). Therefore, the equations solved on grid (k-1) are 

L^^-’+zr^-Ly) 

+(l*-7‘->‘-f‘- 1 ) <2I) 

where /* _1 is the restriction operator. Alternatively, l\_ x is 
the prolongation operator. 0 k is then updated as 

<t>L = <t>L + f L(<P 1 -' - K"<t>L) (22) 

Note that only the change from the previous value 
(<p k ~' -It'^ou) is prolongated to grid k and not the value 

<p k 1 itself. The advantage of using FAS over the Correction 
Scheme is that the solution vector from the fine grid, and not 
just the residuals are transferred to the coarser grids. If 
multiple iterations are performed on a coarse grid, the 
nonlinear operator and the source terms are continuously 
updated. 

The grid iterations can be arranged in a variety of ways 
which affect the overall rate of convergence. The fixed cycle 
is preferred here over an adaptive cycling strategy since it is 
not always possible to assign an optimal smoothing rate as is 
required in an adaptive strategy. The present study employs 
a fixed (1,1) V-cycle with one iteration performed on u, v, 
w, and P per grid during the downward leg and one iteration 
during the upward leg. 


*p = a s + + a* + a s + a H + a L 
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After the pressure correction equation is solved, the cell face 
Cartesian velocities are corrected along with the cell centered 
pressure. The mass fluxes at the cell faces are then 
calculated from the corrected Cartesian velocities. 

M ultigrid Acceleration 

The segregated procedure described has been 
incorporated within a multi grid acceleration technique. A 
full multigrid-full approximation scheme (FMG-FAS) using 
a (1,1) V-cycle has been used. Various multigrid schemes 
and cycles have been proposed and tested in the literature 
(Brandt, 1977; Ghia et ah, 1982; Hortmann et al., 1990, 
etc.). The correction scheme (CS) is the most simple for 
linear problems from the conceptual and algorithmic 
standpoint but is less efficient when used for nonlinear 
problems. The FAS scheme has therefore been utilized The 
full multigrid procedure involves obtaining an approximate 
solution to the problem on coarse grids and subsequently 
prolongating the approximation to finer grids so as to 
provide a good initial estimate. 

The equation set on the finest grid (k) can be written as 
L k <t> k =F k (20) 

Here, L is the nonlinear operator consisting of convection 
and diffusion terms, 0 is the solution vector, and F 
represents the source terms. In the FAS procedure, the 
values calculated on a coarser grid (k-1) are not simple 
corrections to the values on grid (k), instead they are 


Consistent transfer of residuals and solutions is a major 
aspect, of any multi gnd-based algorithm. In the present 
algorithm, cell-center quantities have been restricted using 
full weighting of neighbor values. Restriction of the 
momentum residuals has been performed via summation. 
Corrections have been prolongated by use of trilinear 
interpolation. These restriction and prolongation operators 
are much more straightforward than those used in a 
collocated procedure because of the absence of the 
momentum interpolation. The operators arc essentially the 
same as those for a single component staggered mesh 
procedure (Vanka, 1986a). 

IV. Results 

In the present study, three problems have been 
considered in assessing the performance of the numerical 
method: laminar flow in a curved pipe, laminar flow in a 
combustion chamber, and turbulent flow in a combustion 
chamber. The simulations have been performed with 
multiple velocity grid staggering and with traditional grid 
staggering for purposes of comparison. The results obtained 
using the flux limited quadratic interpolation are also 
discussed. 

Laminar flow in a curved pipe 

The laminar flow in a curved pipe has been considered 
because of the pressure-velocity decoupling which occurs in 
the downstream section of the pipe when conventional grid 
staggering is used. The flow at a Reynolds number of 500 
and radius ratio of 0.217 has been simulated. The inlet plane 
was specified as fully developed pipe flow. The geometry 
used in the simulations is given in Figure 3. A downstream 
pipe section of two diameters was used so as to ensure that 
no flow separation would be present at the fluid outlet 

The simulations performed using the multiple velocity 
grid staggering were observed to yield a grid independent 
rate of convergence (Figure 4). The solution was obtained 
to a normalized mass residual of 0.01% in 
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approximately 50 fine grid iterations. It is apparent from the 
convergence depicted in Figure 5 that the single grid method 
requires many more iterations for convergence than does the 
multigrid method. 



Figure 3: Schematic of curved pipe geometry. 
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Figure 4: Multigrid convergence for calculations of laminar 
flow in a curved pipe. 
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Figure 5: Comparison of muitigrid and single grid rates of 
convergence for curved pipe calculations. 


The single grid convergence histories for the calculations 
using traditional grid staggering are given in Figure 6 for 
different grids. It can be observed that the rate of 
convergence is inferior to that associated with use of multiple 
velocity grid staggering. Table 1 gives the CPU times and 
number of iterations required in obtaining converged 
solutions. All CPU times are for the Cray YMP at the 
NASA Lewis Research Center. It is apparent that although 
the multiple velocity staggering involves the solution of more 
quantities, its CPU time is less than that associated with 
traditional grid staggering. This is because the multiple 
velocity grid staggering is not plagued by pressure-velocity 
decoupling in the downstream section of the pipe. 
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Figure 6: Single grid convergence histories when traditional 
grid staggering is used in the calculation of flow in a curved 
pipe. 

The solution of the curved pipe flow has been compared 
with the experimental results of Enayet et al. (1982). In the 
tests reported in this section, hybrid differencing has been 
used. It can be seen from Figures 7 and 8 that a small 
difference does exist between experiment and predictions. 
The solutions do however indicate that the numerical method 
correctly predicts the flow in this benchmark problem. 
Improved solutions may be obtained by further grid 
resolution or the incorporation of higher order differencing. 



0 0.4 0.8 1.2 1.6 : 

Figure 7: Stream wise velocity profile in the symmetry plane 
at 0=30 degrees. 
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Grid Staggering ■ ||||f 


Iterations CPU time (s) 

Iterations CPU time (s) 

Multi grid 

47 507.4 


Single Grid 

507 2277.4 

1073 3040.2 

I J InminUr flrtll/ ifl ! 


Table 1: Convergence 
curved pipe at Re=500. 
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Figure 8: Stream wise velocity profile in the symmetry plane 

at 0=60 degrees. 

Laminar flow in the dilution zone of a gas turbine 
combustion chamber 

Laminar flow in the dilution zone of a gas turbine 
combustion chamber was calculated as a test of the numerical 
method. For laminar flow, two cases have been considered: 
isothermal flow and laminar mixing of hot and cold flowing 
streams. The geometry considered is the same as that used 
by McGuirk and Palma (1992) and is depicted in Figure 9. 
For the isothermal flow, the inlet Reynolds number was 
equal to 27 with the top jet issuing into the chamber at a 
Reynolds number of 39 and the bottom jet issuing at a 
Reynolds number of 23. The laminar mixing problem 
considered an inlet Reynolds number of 27, a top jet at 
Reynolds number of 117, and a bottom jet at Reynolds 
number of 69. The inlet temperature for the mixing problem 
was 900 K while both jets were maintained at a uniform 
temperature of 300 K. The isothermal flow will be 



Figure 9: Schematic of the dilution zone of a gas turbine 
combustion chamber. 


From Figure 10 it is apparent that the numerical method 
exhibits a grid independent rate of convergence.Asolunon 
to 0.1% normalized mass residual is obtained in 
approximately 50 fine grid iterauons. c °mpanson of the 
multiple velocity scheme with the traditional grid staggenng 
scheme reveals that both methods converge at the same rate 
for both the single grid (Figures ; 11 and 12) andthemulu^d 
formulations (Figures 10 and 13). However, *«£PUnnres 
required for the two methods are different. From Tables 2 
and 3 it is seen that the time required by traditional grid 
staggering is approximately 2/3 of the mulup e velocity 
scheme. Thus, in order for the multiple velocity gnd 
staggering to be beneficial to this problem, it must converge 
in fewer than 2/3 the number of iterauons required by the 
traditional grid staggering scheme. This is not seen to be the 
case. Tables 2 and 3 indicate that when using the muluple 
velocity grid staggering or the traditional grid staggering, the 
muldgrid convergence is faster than the corresponding single 
grid convergence by a factor of 2.5. 



Fieure 10: Multigrid convergence for laminar, isothermal 
flow in a combustor when multiple velocity gnd staggenng 
is used. 

10° 



Figure 11: Single grid convergence for laminar, isothermal 
flow in a combustor when multiple velocity grid staggering 
is used. 
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Figure 12: Single grid convergence for laminar, isothermal 
flow in a combustor when traditional grid staggering is used. 
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Figure 13: Multigrid convergence for laminar, isothermal 
flow in a combustor when traditional grid staggering is used. 

The mixing problem also reveals that both grid 
staggering schemes have the same rate of convergence 
(Figure 14 and 15). Therefore, the traditional grid 
staggering requires less work as seen in Table 2. The CPU 
time required by the multiple velocity grid staggering scheme 
is larger by a factor of 1.6. In the muldgrid mode, a strict 
grid independent rate of convergence was not obtained as 
evidenced by Figures 16 and 17. The observed behavior is 
due to the fact that the energy equation has not been 
incorporated into the V-cycle and thus its convergence is not 
accelerated. In order to obtain a striedy grid independent rate 
of convergence all equations should be contained in the V 
cvcle. 

10 l r 
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Figure 14: Single grid convergence for laminar mixing flow 
in a combustor when multiple velocity grid staggering is 
used. 
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Figure 15: Single grid convergence for laminar mixing flow 
in a combustor when traditional grid staggering is used. 
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Figure 16: Multigrid convergence for laminar mixing flow 
in a combustor when multiple velocity grid staggering is 
used. 
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Figure 17: Multigrid convergence for laminar mixing flow 
in a combustor when traditional grid staggering is used. 
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51 

129.4 

54 | $2.5 1 

| Cold mixing 

144 

403.5 



Table 2: Multigrid convergence characteristics for laminar flow in the dilution zone of a gas 
turbine combustion chamber. 



| Multiple Velocity Grid Staggering 



■Hsnssm 

CPU time (s) 


BES1SES3B1B 

Isothermal 

329 

3431 “ 

295 

177.9 

Cold mixing 

522 

§7071 ~ i 

523 

405.2 


Table 3: Single grid convergence characteristics for laminar flow in the dilution zone of a 
gas turbine combustion chamber. 


The solutions obtained are shown in Figures 1 8-22. It is 
apparent that the opposed jets merge for the case of a cold jet 
mixed with a hot fluid. The isothermal flow does not exhibit 
this jet merging. This is to be expected since the cold jets 
which issue into the chamber in the mixing flow case have a 
much larger momentum. The maximum momentum ratio for 
the isothermal case is 6.9, whereas the maximum momentum 
ratio in the mixing flow problem is 20.8. Also, in the 
mixing flow case a forward recirculation zone forms in the 
region where the opposed jets merge. This is not observed 
in the case of isothermal flow where the jets do not merge. 
Selected temperature distributions at various planes in the 
combustor are shown in Figure 20. The temperature 
distribution in the symmetry plane (Figure 20a) indicates that 
the high temperature flow from the inlet experiences more 
cooling along the bottom boundary than along the top. 
Figure 20b reveals that the fluid in the injection plane is 
advected toward the upper wall. These observations are 
further verified by Figure 20c which shows the temperature 
distribution in the exit plane of the combustor. The fluid at 
the top wall of the outflow nozzle is at a higher temperature 
than at the bottom walL 



Figure 18: Velocity field in injection plane for laminar, 
isothermal flow in a combustor. 



Figure 19: Velocity field in the symmetry plane for laminar, 
isothermal flow in a combustor. 



Figure 21: Velocity field in injection plane for laminar 
mixing flow in a combustor. 



Figure 22: Velocity field in the symmetry plane for laminar 
mixing flow in a combustor. 


Turbulent flow in the dilution zone of a gas turbine 
combustion chamber 

The third case considered is the isothermal turbulent flow 
in the dilution zone of a gas turbine combustion chamber for 
an inlet Reynolds number of 12,000. The maximum jet-to- 
inlet momentum ratio is 6.9. The multigrid solution of this 
problem indicates that some deterioration in the rate of 
convergence occurs with grid refinement. As in the case of 

the laminar mixing, the k and e equations have not been 
incorporated into the multigrid V -cycle. Several researchers 
have dealt with the issue of incorporating the k-e turbulence 
model into a multigrid solution procedure (Vanka, 1987; 
Rubini et aL, 1992; Shyy et al., 1993), but the rates of 
convergence for turbulent flow have still not been shown to 
be equivalent to that which is observed for laminar flows. 
The solution to 0.2% normalized mass residual was 
accomplished in approximately 150 fine grid iterations using 
both the multiple velocity grid staggering scheme as well as 
traditional grid staggering (Figures 23 and 24). The single 
grid convergence was obtained in approximately 555 
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CPU time (s) 

■BTTTTTifrTrHl 

CPU time (s) 


150 


150 

383.2 


559 

1214.3 | 

555 

930.3 


Table 4: Convergence times and total number of iterations required for calculation of 
isothermal turbulent flow in the dilution zone of a gas turbine combustor. 


iterations using both grid staggering schemes. Here also the 
traditional grid staggering requires a smaller CPU time (see 
Table 4). 

The solution of the isothermal turbulent flow reveals that 
the opposed jets do impinge for this case. A small 
recirculation zone forms in front of the impingement region 
just as in the laminar mixing case (Figure 25). The 
recirculation region aft of the bottom jet which was present 
for laminar flow (Figure 18) is no longer present for 
turbulent flow. The flow in the symmetry plane between 
combustion chambers (Figure 26) indicates that the fluid 
from either side of the chamber impinges at this plane and 
advects outward toward the walls. 



Figure 23: Multigrid convergence for calculation of 

turbulent flow in a combustor when multiple velocity grid 
staggering is used. 



Number of iterations 

Figure 24: Multigrid convergence for calculation of 

turbulent flow in a combustor when traditional grid 
staggering is used. 



Figure 25: Velocity field in injection plane for turbulent flow 
in a combustor. 



Figure 26: Velocity field in the symmetry plane for turbulent 
flow in a combustor. 


Results obtained with flux limited quadratic interpolation 

The laminar flow in a curved pipe of Reynolds number 
500 and radius ratio 0.217 has also been calculated with a 
flux limited form of QUICK. The calculations have been 
made with the multiple velocity grid staggering scheme. The 
convergence history on a 60x32x16 grid using three levels in 
a (1,1) V-cycle is shown in Figure 27. The convergence is 
compared with that of the first order upwind scheme. It is 
observed that the higher order discretization converges only 
to 0. 1 % after which a high frequency oscillation is exhibited. 
This behavior was also seen in the single grid mode as well 
as when the multigrid cycling was switched off at 1% 
residual. A similar behavior was also observed for flow in a 
cube with a moving top wall. For this case, the limiting 
residual decreased as finer grids were considered. For a 
10x10x10 grid the limit residual was approximately 0.1% 
while for a 40x40x40 grid, the method converged to 0.01% 
normalized mass residual. The flux limited quadratic 
interpolation is a nonlinear discretization which changes with 
the local solution. It is speculated that the limit residual may 
be the result of oscillation in the deferred correction source 
term arising from the flux limiting process. 
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Figure 27: Convergence history of laminar pipe flow 
calculation on a 60x32x16 grid using two different 
discretizations. 


V. Conclusions 

A multigrid method for internal flows in complex 
geometries based on multiple velocity grid staggering has 
been described. The numerical method has been tested for 
laminar flow in a curved pipe and for laminar and turbulent 
flow in the dilution zone of a gas turbine combustion 
chamber. The convergence characteristics for the method 
have been compared to the convergence obtained for use of 
traditional grid staggering. It has been observed that the use 
of the multiple velocity grid staggering arrangement is 
superior to traditional grid staggering only when a ninety 
degree bend exists in the problem geometry. In this case, 
the multiple velocity arrangement avoids the pressure- 
velocity decoupling phenomena which occurs when 
traditional grid staggering is used for such problems. For 
the simulation of flow in the dilution zone of a gas turbine 
combustion chamber, it was observed that both grid 
staggering arrangements required approximately the same 
number of iterations for convergence. In this case, the 
traditional grid staggering is to be favored due to its lower 
work count per iteration. One potential advantage of the 
multiple velocity grid staggering arrangement involves the 
use of coupled solvers for the solution of the momentum and 
continuity equations. The multiple velocity grid staggering 
arrangement facilitates use of a symmetrical coupled Gauss- 
Seidel operator (Vanka, 1986c). It is not feasible to use 
such a coupled solver if only one Cartesian velocity 
component is stored at each cell face. The potential benefit 
arises from the lower work count of the SCGS operator in 
comparison to the SIMPLE algorithm used in this study. 
The SIMPLE algorithm requires approximately 30% more 
work than one sweep of SCGS (Sockol, 1993). If it is 
possible to obtain a converged solution with the SCGS 
operator in the same number of iterations as the SIMPLE 
algorithm when a multiple velocity staggered grid is used, 
the cost incurred by having more Cartesian velocity 
components may be alleviated. 
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Figure 20: Temperature distributions for laminar mixing 
flow in a combustion chamber: (a) lengthwise plane at 0=0 
degrees (b) cross- stream plane at injectors (c) exit plane. 
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ABSTRACT 

The present study compares coupled and segregated 
relaxation operators with multigrid acceleration for the 
calculation of jet impingement now fields. The symmetrical 
coupled Gauss-Scidcl (SCGS) scheme and the SIMPLE 
algorithm have been used as relaxation operators. The coupled 
operator is observed to require less work for convergence than 
the segregated method. When additional scalar equations were 
solved, the segregated method exhibited a definite 
deterioration in convergence rate. The coupled method 
revealed similar behavior as regards the solution of turbulence 
equations. Unlike the segregated operator, the convergence of 
the coupled operator did not deteriorate with the addition of the 
energy equation for high temperature flow calculations. The 
calculations performed in these tests required cell aspect ratios 
as large as sixteen in order to resolve important flow features. 
The SCGS operator provided rapid convergence even when 
these aspect ratios were used. 

INTRODUCTION 

In recent years, the multigrid technique has been shown to 
significantly improve the performance of iterative 
computational algorithms for the solution of the discretized 
fluid flow equations (Brandt, 1980; Demurcn, 1989; Vanka, 
1986; Shyy el al., 1993). In the multigrid technique, a number 
of coarse grids are used in conjunction with a given fine grid in 
order to accelerate the rate of convergence of the iterative 
procedure. Residuals lhai arc slow io converge on a given grid 
arc interpolated and solved on progressively coarser grids. 
The corrections implied by the solutions on the coarser grids 
are then applied to the fine grid solution. The primary 
advantage of the multigrid method is that the rate of 
convergence can be made independent of the number of discrete 


nodes in the computational domain and the work count can be 
made to vary as 0(n). In contrast, the convergence of single 
grid algorithms deteriorates as the number of grid nodes is 
increased. The cause of such a behavior is attributed to the 
slow convergence of the low frequency errors that are present 
in any partially converged solution. By restricting such errors 
to coarser grids where they appear as high frequency errors, 
their rate of convergence can be improved. Beginning with 
the pioneering works of Brandt (1977), the multigrid method 
has been demonstrated to work well in a variety of 
applications including the computation of viscous internal 
fluid flows. 

In any multigrid scheme, the relaxation procedure used to 
resolve the errors on a given grid plays an important role in 
the success of the overall numerical method. For viscous 
internal fluid flows, several relaxation procedures can be 
designed. Several computational parameters play an important 
role in the relative behavior of the relaxation scheme. These 
include the geometrical configuration, the aspect ratio of the 
grid cells, the speed of the flow (supersonic vs. subsonic), the 
presence of interior flow obstructions and the flow Reynolds 
number. It is not possible to provide a priori information on 
the best relaxation procedure for a given flow because of the 
complex nonlinear nature of the equations. As a result, much 
of the information must be gained through the testing of their 
performance in actual implementations. 

For viscous internal flows dominated by the strong 
coupling between the pressure and the velocity fields, two 
main relaxation strategics have been advocated. In one of 
these, the momentum and continuity equations are solved in a 
coupled manner, thus preserving the strong pressure-velocity 
coupling. In the other, the momenium and continuity 
equations are solved in a segregated manner, wherein the 
momentum equations arc first solved with an approximate 
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pressure field. The pressure Field is then corrected using the 
local imbalances in the mass continuity equation. A pressure* 
correction equation is derived by substituting approximate 
forms of the discrete momentum equations into the continuity 
equation. Several pressure update strategies have been 
proposed in the past, having built on the original concept of 
the SIMPLE algorithm of PatanJkar and Spalding (1972). 

Recent progress in the development of coupled operators 
has resulted in several procedures which have been observed to 
be more computationally efficient than segregated operators 
on certain problems. Galpin et al. (1985) developed the 
coupled equation line solver (CELS) which simultaneously 
solves for the velocities and pressures along a given grid line. 
In tests with two improved versions of the SIMPLE algorithm 
(SIMPLEC and SIMPLER), it was observed that CELS was not 
as sensitive to the choice of relaxation factor as was SIMPLEC 
and SIMPLER. The observation was a significant one since 
the optimal relaxation factor is not known a priori. The 
symmetrical coupled Gauss-Seidcl (SCGS) operator proposed 
by Vanka (1986) is a point solver in that velocities and 
pressures for only one cell are coupled. One advantage which 
is observed for point solvers as compared to line solvers is a 
lower work count per iteration. Dcmuren (1989) extended the 
SCGS operator to solve implicitly all the pressures along a 
line (CLSOR) or a plane (SIM) with subsequent update of the 
velocities at each cell. The SCGS, coupled line successive 
over-relaxation (CLSOR), and strongly implicit (SIM) 
relaxation operators all involve local coupling between 
velocity and pressure. The extension of SCGS to CLSOR and 
SIM provided additional coupling for pressures along a line or 
plane. In comparisons of SCGS, CLSOR, SIM, and SIMPLE 
for a lid driven cavity flow, Dcmuren found that the SCGS 
operator gave the best performance of all methods in terms of 
CPU time and computer memory when implemented in the 
context of a multigrid method. The observation was made for 
nearly isotropic coefficients in the difference equations. The 
use of large cell aspect ratios (Ax/Ay) leads to anisotropy in 
the coefficients of the discrete equations. It was observed that 
SCGS did not yield the best performance when large aspect 
ratios were used. A similar observation was made by Rodi et 
al. (1987) as regards simulation of three-dimensional 
problems involving large aspect ratios. Recently, Sockol 
(1993) has compared SCGS, CLSOR, fully coupled line block 
Gauss Seidel (CLBGS), and SIMPLE relaxation procedures in 
the context of a mulligrid method. The comparisons weTC made 
for three different flows: lid driven cavity flow, developing 

channel flow, and open cavity flow. In tests at various 
Reynolds number and cell aspect ratio, it was concluded that 
SCGS offered the best mix of robustness and computational 
speed. It was also observed that SCGS was less sensitive to 
the presence of comer singularities than were the other 
relaxation operators, but that the option of using CLBGS 
should be retained for regions of strongly aligned flow (e.g. 
flow in a straight pipe). 

The present work is motivated by the desire to efficiently 
solve the flow fields created by the exhaust jets of a Short Take 
off and Vertical Landing (STOVL) Aircraft hovering in ground 


proximity. The determination of these flowfields is important 
to the estimation of the temperature of the hot gases ingested 
by the engine and to the determination of the loss of lift caused 
by vortices underneath the aircraft. Such flows provide several 
problems for any numerical method. Since jet impingement 
phenomena are inherently three-dimensional, a large number 
of cells must be used to coveT the flow domain. As noted, 
convergence deteriorates with gnd resolution for a single grid 
method, hence a multi-grid method is extremely important for 
such three-dimensional calculations. The use of fine grids in 
the jet impingement region to resolve sharp gradients, gives 
rise to large aspect ratios in the region of the ffeestream where 
no such resolution is needed. An additional problem, common 
to any turbulent flow, is that of equation coupling. The 
coupled methods which have been discussed involve coupling 
between velocities and pressures. The solution of scalar 
equations such as k and C is usually accomplished in a 
segregated manner. 

In the present study, a comparison of the segregated and 
coupled strategies has been made for the computation of the 
flow fields created by twin impinging jets in the presence of a 
head wind. The performance of ihe two relaxation operators is 
evaluated in terms of convergence and CPU time. A 
description of the governing equations and discretization is 
given first. The segregated and coupled procedures are then 
described followed by a discussion of the multigrid 
implementation of the two relaxation operators. The results 
obtained from the numerical experiments indicate that the 
coupled operator used in this study (SCGS) is computationally 
more efficient than the segregated procedure (SIMPLE) even 
when large cell aspect ratios are present 

GOVERNING EQUATIONS AND 

DISCRETIZATION 

The continuity, momentum, energy, and turbulence 
equations may be conveniently written in terms of a general 
conservation equation which represents a balance of 
convection, diffusion, and source terms (Patankar, 1980): 


— (pu<p) + — (pv<t>) + — (p”4>) = 
dx dy dz 



r* and S ^ are the diffusive exchange coefficient and source 
term respectively, for the general variable 0. The values of 

and S ^ associated with a given flow variable are given in 
Tabic 1. The strong conservation law form of the governing 
equations has been retained so as to facilitate use of a finite 
volume discretization. 

For turbulent flows, the Reynolds averaged equations are 
closed by use of a k-£ turbulence model (Launder and Spalding, 
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1974). The model relates the Reynolds stresses to the mean 
velocity field through use of the Boussinesq approximation. 

The associated eddy viscosity is given by p.j = pC^k^ Je 

where modeled equations for k and £ are solved. The turbulence 
model used in the present study incorporates the wall functions 
of Launder and Spalding (1974) to specify boundary conditions 
for the flow variables near a solid boundary. 

The Finite volume method is used to discretize the 
governing equations. In this method, the integral and not the 
differential form of the governing equations is approximated. 

In this method, the computational domain is divided into 
numerous control volumes. The integral form of the 
conservation laws is obtained for each control volume by 
integrating the governing equations (Eqn. 1) over the bounds 
of each control volume. The velocity and pressure cells have 
been staggered relative to one another so as to avoid the 
pressure-velocity split which occurs for simulations of 
incompressible flows (Harlow and Welch, 1965). The discrete 
expression for the combined convective/diffusive flux utilizes 
the hybrid differencing procedure of Spalding (1972). Having 
discretized the governing equations, the problem becomes one 
of solving for a set of coupled, nonlinear algebraic equations. 
The nonlinear algebraic equations are linearized by the 
procedure of lagged coefficients (Anderson et al., 1984) and 
solved through the process of iteration. The discrete equations 
are written in the following form: 

A P*P = 1 A nb*nb + Az (2) 

In order to maintain numerical stability during the 
iterative process, it is necessary to dampen the successive 
changes of the flow variables. Under-relaxation has been 
implemented implicitly by changing the discrete equations as 
follows (Patankar, 1980): 


= lA nb/ a 

(3) 

M 

l 

a 

5: 

(4) 


where a is the relaxation factor. The optimal relaxation factor 
varies with the problem solved and the characteristics of the 
finite -difference grid. 

SOLUTION METHODOLOGY 

The present study involves a comparison of the 
performance of multigrid based segregated and coupled 
solution procedures for the simulation of jet impingement 
flows. The segregated approach used in this investigation is 
the SIMPLE algorithm of Patankar and Spalding (1972). The 
coupled solution procedure used is the symmetrical coupled 
Gauss-Seidel (SCGS). proposed by Vanka (1986). The 
segregated and coupled solution procedures arc embedded in a 
multigrid cycle which serves to accelerate their rates of 
convergence. 


Segregated Method 

In the segregated relaxation method, each equation is 
solved in a sequential fashion. On a given grid, the following 
procedure is used to iteratively solve for the flow variables: 

1 ) Solve for the u component of velocity 

2) Solve for the v component of velocity 

3) Solve for the w component of velocity 

4) Solve for the pressure correction. F 

5 ) Update u, v, w, and P based on the solution for F 

6) If on the locally finest grid, solve for scalar variables 
(T,k,e) 

where the pressure correction equation is written in a form 
analogous to the velocity and scalar equations: 

A p P 'p=' LA nb P 'nb +Sm ( 5 ) 

The quantity S m represents the mass imbalance for an 
individual cell. The iterative process for any particular flow 
variable is not carried to convergence since subsequent 
iterations on the other flow variables will change the 
coefficients of the discrete equation. In general, five iterations 
of an alternating line Gauss-Seidel procedure were performed 
for each velocity and scalar variable while fifteen iterations 
were performed on the pressure correction equation. The six 
steps given above constitute one "outer" iteration, i.e. 
iterative update of all flow variables. A numerical method 
employing only a single grid would simply continue these 
outer iterations until the residuals in the governing equations 
were reduced below a prescribed level. In the present multigrid 
method, the outer iterations are performed on a given grid prior 
to restriction or prolongation to another grid. It is important 
to note that during the execution of a V-cycle, the scalar 
equations are only solved on the locally finest grid. The 
values of T, k, and £ on the finest grid are held fixed on the 
coarse grids during the V-cycle. The details pertaining to the 
multigrid acceleration technique are described below. 

The segregated solution procedure has been shown to be 
convergent for a wide range of flow problems although it can 
be very slow to converge on a single grid (Rayner, 1991). One 
potential shortcoming of the segregated method is the 
pressure-velocity coupling. The segregated method involves 
writing the continuity equation as a pressure correction 
equation. There are various means of deriving a pressure 
correction equation. The different pressure correction 
equations do not affect the final solution since all pressure 
corrections will be zero when the solution has converged. The 
rate of convergence however can be different for the various 
pressure correction practices. The coupled solution procedure 
seeks to improve on the coupling between pressure and 
velocity and it eliminates any ambiguity involving pressure 
correction equations simply because the primitive variable 
form of the continuity equation is retained. 
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Coupled Method 

A coupled solution procedure solves for the velocities end 
pressures simultaneously. This is in contrast to the segregated 
procedure which solves these quantities in a sequential manner. 
The coupled solution seeks to attain a better coupling between 
velocities and pressures by retaining the primitive variable 
form of the continuity equation. A pressure correction 
equation is not needed in a coupled solution method. The 
coupled method used in the present investigation is the 
symmetrical coupled Gauss-Seidel (SCGS) proposed by Vanka 
(1986). The method was chosen because of the strong local 
pressure-velocity coupling as well as for its low work count. 
The SCGS operator also requires less storage than the coupled 
line or plane solvers. This can be a substantial benefit when 
very fine grids are required in 3-D calculations. 


The SCGS operator considers a local coupling between 
pressure and velocitie*. For a 3-D geometry, the six 
momentum equations for a given finite volume, along with the 
continuity equation provide the seven equations necessary to 
find the six components of velocity located at the cell faces 
and the cell centered pressure. These equations are written in 
terms of current approximations for the unknowns plus a 
correction: 

old 

u = u + U 


old 

V = V + V 

old 

w = w + w 


( 6 ) 


p = p old + P' 


The following 7x7 matrix can be analytically inverted to 
obtain the corrections to the local unknowns: 


( A P 1-1/2 

0 

0 

0 

0 

0 

SySz 


“«-V 2 


IV 1 

0 

( A p)i+\/2 

0 

0 

0 

0 

-SySz 
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0 

0 

wu 

0 

0 

0 

SxSz 
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0 

0 
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0 

0 
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X 
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0 

0 
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0 

0 

V ^'> + 1/2 
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wu 

0 

0 

wu 

Sxdy 
- SxSy 


w 'k- 1/2 
W k+\J7 

J 


R k- 1/2 
R k+ 1/2 
_ . 

-SySz 

SySz 

- SxSz 

6x5? 

- 5xSy 

SxSy 

0 



where is the central coefficient in the discrete equation for 
the variable 0 and is the residual (Vanka, 1986). The seven 
variables at each cell are solved for using algebraic relations. 
The computational domain is swept through in a lexicographic 
manner in order to update all velocities and pressures. When 
the entire domain has been swept through, each velocity will 
have been updated twice and each pressure once. 


Multlarld Acceleration 

The segregated and coupled solution procedures described 
have been incorporated within a multigrid acceleration 
technique. A full approximation scheme-full multigrid (FAS- 
FMG) algorithm using a (1.1) V -cycle has been used. Various 
multigrid schemes and cycles have been proposed and tested in 
the literature (Brandt, 1977; Ghia et al., 1982; Hortmann et 
al., 1990, etc.). The FAS method has been used since the 
problem is nonlinear. The correction scheme (CS) is 
conceptually and algorithmically more simple for linear 
problems, but is less efficient when used for nonlinear 
problems. The full multigrid procedure involves 
approximately solving the problem on coarse grids and 


r~ 

subsequently prolongating the solution to the finer grids so as 
to provide a good initial estimate. 

The equation set on the finest grid (k) can be written as 

Z.V = F k (8) 


Here, L is the nonlinear operator consisting of convection and 
diffusion terms, 0 is the solution vector, and F represents the 
source terms. In the FAS procedure, the values calculated on a 
coarser grid (k-1) are not simple corrections to the values on 
grid (k); instead they are approximations on grid (k-1) to the 
correct values on grid (k). Therefore, the equations solved on 
grid (k-1) are 


1 = F k ~ } + l k Vf* - f-V) 

v (9) 

k-1 k 

where is the restriction^ operator. Alternately, ^ is 

the prolongation operator. 0 is then updated as 


* k 

“new 


*old 




k-l 




(10) 
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Note that only the change from the previous value 
is prolongated to grid k and not the value 

itself. The advantage of using FAS over the Correction 
Scheme is that the solution vector from the fine grid, and not 
just the residuals are transferred to the coarser grids. 
Additionally, if multiple iterations are performed on a coarse 
grid, the nonlinear operator and the source terms are 
continuously updated. 

The grid iterations can be arranged in a variety of ways 
which affect the overall rate of convergence. The fixed cycle is 
preferred here over an adaptive cycling strategy since it is not 
always possible to assign an optimal smoothing rate as is 
required in an adaptive strategy. The present study employs a 
fixed (1.1) V-cycle with one iteration performed on u, v, w. 
and P per grid during the downward leg and one iteration during 
the upward leg. 

Consistent transfer of residuals and solutions is a major 
aspect of any mulligrid-based algorithm. In the present 
algorithm, cell-center quantities are restricted using full 
weighting of neighbor values. Restriction of residuals is 
accomplished by requiring that the per unit volume residual on 
the fine grid is maintained on the coarse grid. Cell-face 
velocities are restricted so that mass conservation is 
maintained on both the fine and coarse grids. Corrections are 
prolongated by trilinear interpolation. 


RESULTS 

The simulation of jet impingement flows using coupled 
and segregated solution procedures has been accomplished for 
four test cases: laminar cold and hot impinging jets, turbulent 
cold and hot impinging jets. The basic flow fields observed 
for the turbulent flow cases will be briefly described as these 
are more relevant to the STOVL problem than the laminar flow 
cases. The convergence characteristics of the two methods are 
of primary importance in this investigation. The number of 
iterations required for convergence as well as the overall CPU 
lime required for solution on a Cray Y-MP are considered in the 
analysis of the two methods. 


Problem Definition 

The geometry considered for all test cases is depicted in 
Figure 1. Two jets impinge normal to an infinite flat plate. 
The parallel plate is at a distance of five jet diameters from the 
jet exhaust and the two jets are placed six diameters apart. A 
cross flow is specified in a direction normal to the jets and has 
a magnitude of 3% relative to the jet velocity. Due to 
symmetry of the flow field in the z direction, only one half of 
the jet flow field is actually simulated. The computational 
domain was 69 jet diameters in the streamwise direction (x), 5 
diameters in the cross stream direction (y), and 12.5 diameters 
in the spanwisc direction (z). These dimensions were deemed 
necessary in order to accurately represent the flow 
characteristics present in an impinging jet flow without 
having interference from the simulated boundaries. The largest 


cell aspect ratio was sixteen for use of a fine grid of 
120x20x36. Other details related to the geometry and flow 
variables used are given in Tables 2 and 3. 

The simulation of any flow field requires that appropriate 
boundary conditions be specified. At the inlet of the domain, 
uniform velocity and temperature distributions have been 
specified. No slip, no penetration adiabatic walls have been 
specified at the solid surfaces. The outflow boundary was 
treated as fully developed flow in the present simulation, i.e. 


3u 

dv 

d* 

dh 

— 

= = 

= 

■" 

dx 

dx 

dx 

dx 


(11) 


The z-minus boundary bisected the exhaust jets and therefore a 
symmetry boundary condition was specified. The z-plus 
boundary has been treated as a solid surface in the present 
calculations. The jets themselves were represented as mass, 
momentum, and energy sources and the profiles of velocity, 
density, temperature, k, and £ were specified as uniform. 


Convergence and Execution TinfeS 

The simulations of impinging jet flow have been 
performed on a finest grid of 120x20x36 cells with a total of 
thjee (3) grids used in the solution. The relaxation factors and 
number of iteration sweeps used on the scalar equations for the 
four test cases are given in Table 4. 

The convergence histories for the laminar, cold jet 
simulation using coupled and segregated operators are given in 
Figures 2 and 3. In Figure 2, the convergence history of the 
simulation as performed on all three grids using SCGS is 
presented. A similar convergence history is given in Figure 3 
for the SIMPLE algorithm. It is observed from these figures 
that the multigrid rate of convergence is obtained with both 
relaxation operators. In other words, as the grid is refined 
from 30x5x9 out to 120x20x36, the convergence rate does not 
deteriorate. The number of iterations required for convergence 
on the coarse grid is essentially the same as the number of 
iterations required for convergence on the finest grid. This 
grid independent convergence is not achieved with single grid 
procedures and is the primary motivation for use of multigrid 
methods. It is also apparent that the coupled relaxation 
operator requires approximately 10-20 fewer iterations for 
convergence than does the segregated relaxation operator. 

Table 5 reveals that the actual CPU time required for 
convergence is smaller for the coupled operator. This fact is 
not only the result of fewer iterations, but also due to the lower 
work count for the SCGS operator when compared to the 
SIMPLE method. The SCGS operator requires less 
computational work per iteration than SIMPLE. Sockol 
(1993) reports that SIMPLE requires approximately 30% more 
work per iteration than SCGS. 

Comparison of convergence rates for the laminar and 
turbulent flow problems reveals that almost twice as many 
iterations (and hence work) is required for a turbulence 
simulation. This is understood when one considers that the 
scalar equations are not incorporated into the V-cycle of the 


98 



multigrid procedure. In particular, it is the turbulence 
equations which slow the rate of convergence since in the 
laminar, hot jet problem the energy equation is solved in 
single grid mode without a significant deterioration in the rate 
of convergence for use with the coupled operator (compare 
Figures 2 and 4). When the segregated relaxation procedure is 
used, solving the energy equation in single grid mode is 
sufficient to cause a deterioration in the rate of convergence 
(see Figures 3 and 4). These observations are also made for 
turbulent flow, where the solution of scalar equations seems to 
be more deterimental to the convergence rate when a 
segregated operator is used. Note from Table 5 that turbulent 
hot jet simulations performed with the segregated method 
required two more CPU minutes than did the turbulent cold jet 
simulations. In contrast, the turbulent hot jet simulations 
performed with the coupled operator actually required less time 
than the turbulent cold jet simulations. 

From Figures 4 and 6 it is observed that the hot jet 
simulations required more work when the segregated relaxation 
operator was used. In this study, a work unit is defined as 
being an equivalent fine grid iteration. Table 5 reveals that 
almost twice as many work units were required for the 
segregated method as compared to the coupled method for hot 
jet calculations. In this case, not only does the coupled 
operator require less work per iteration, it also requires fewer 
iterations for these calculations. 

Flow Field Description 

The velocity fields in the symmetry and ground planes of 
both cold and hot turbulent flows are given in Figures 7-10. 
The exhaust jets impinge on the ground plane forming wall 
jets which flow radially in all directions from the point of 
impact. A ground vortex forms when the wall jets from the 
fore and aft exhaust jets are stagnated by one another (see 
Figures 7 and 9). It is apparent from Figures 7 and 9 that no 
up wash fountain forms as seen in actual STOVL flow fields 
(Mac Lean, et ah, 1992). In other words, the stagnated wall 
jets do not turn upward and impinge on the underside of the 
exhaust plate. The ground vortex for the cold flowing jet is 
seen to be slightly larger than that of the hot flowing jet. 
Also, the stagnation point of the two wall jets is observed to 
be nearer the fore exhaust jet in the cold flow case. From 
Figure 9, one can observe that for the hot exhaust jet, the 
ground vortex partially blocks the aft exhaust jet. This 
behavior is not observed for the cold flow. The stagnation 
line is bent in the direction of the freestream flow for the hot 
exhaust jet case, while in the cold jet case it is almost 
perpendicular to the symmetry plane. 

CONCLUSIONS 

The present study has compared coupled and segregated 
relaxation operators with multigrid acceleration for the 
calculation of jet impingement flow fields. The symmetrical 
coupled Gauss-Seidel (SCGS) scheme and the SIMPLE 
algorithm have been used as relaxation operators. The coupled 


operator is seen to require less work for convergence than the 
segregated method. When additional scalar equations were 
required to be solved, the segregated method showed a definite 
deterioration in convergence rate. The coupled method 
revealed similar behavior as regards the solution of turbulence 
equations. Unlike the segregated operator however, the 
convergence of the coupled operator did not deteriorate with 
the addition of the energy equation for hot flow calculations. 
In test cases for which both methods required essentially the 
same number of iterations to converge, the lower work count 
per iteration for the SCGS operator led to a smaller execution 
time when compared to use of the SIMPLE algorithm. Finally, 
the calculations performed in these tests required cell aspect 
ratios as large as sixteen in order to resolve important flow 
features. The SCGS operator provided rapid convergence even 
when these aspect ratios were used. Our observations are 
consistent with the finding of Sockol (1993) and indicate that 
the SCGS operator can be used when large aspect ratios are 
required as in the case of impinging jet flows. 
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TABLE 2. GEOMETRICAL AND FLOW PARAMETERS 


Jet exit diameter, D 

3.66 cm 

Jet spacing, S 

21.96 cm 

Domain length, L 

252.54 cm 

Domain height, H 

18.3 cm 

Domain half width, W 

45.75 cm 

Jet exit velocity, V ^ 

100.0 m / s 

Cross stream velocity, V m 

3.0 m / s 


TABLE 3. JET TEMPERATURES AND REYNOLDS NUMBERS FOR THE 
VARIOUS TEST CASES CONSIDERED. 




Tict (K) ] 

Re 

1 

300 

238 

Laminar hot jet 

1000 1 

71 

Turbulent cold jet 

300 1 

238,000 


1000 

70,800 


Laminar 

cold 


Laminar 

hot 


Turbulent 

cold 


TABLE 4. SOLUTION PARAMETERS FOR MOMENTUM AND SCALAR EQUATIONS. 


Segregated Relaxation Procedure 


<*U 

ot p 

dj 

«k,e 

dp 

0.5 

0.3 

NA 

NA 

NA 

0.5 

0.3 

0.9 

NA 

NA 

0.5 

0.3 

NA 

0.5 

0.7 

0.5 

0.3 

0.9 

0.5 

0.5 


Coupled Relaxation Procedure 


a u a p ar a k , e % 


0.5 1.0 NA NA NA NA 


0.5 1.0 0.9 NA NA 3 


0.5 

1.0 

NA 

0.5 

0.5 

3 

0.5 

1.0 

0.9 

0.5 

0.5 

3 
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TABLE 5. COMPARISON OF SEGREGATED AND COUPLED OPERATORS 


Work Unit 


Laminar cold jets 

98.1 

87.8 

271.5 

135.5 

Laminar hot jets 

118.4 

57.4 

362.1 

143.0 

Turbulent cold jets 

242.5 

219.5 

975.3 

719.4 

Turbulent hot jets 

287.1 

150.0 

1108.8 

482.8 








FIGURE 1 . COMPUTATIONAL DO MAW IN THE X-Y PLANE. 
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Grid 2 1 60*10*18 J 
Grid 3 1 120*20*36] 


Number of iterations 

Figure 2 Convergence history for laminar twinjet impingement using 
a coupled relaxation operator. The jet exit temperature was 300 K. 
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Figure 3 Convergence history for laminar twinjet impingement using 
a segregated relaxation operator. The jet exit temperature was 300 K. 
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Figure 4 Comparison of convergence history for laminar twinjet 
impingement The jet exit temperature was 1000 K.. 
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Figure 5 Comparison of convergence history for laminar rwinjet 
impingement The jet exit temperature was 300 K. 
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Figure 6 Comparison of convergence history for turbulent twinjet 
impingement. The jet exit temperature was 1000 K. 
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Figure 7. Velocity distribution in the plane of symmetry (x-y plane) for a turbulent 
impinging jet. The jet temperature is 300 K. 
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Figure 8. Velocity distribution in the ground plane (x-z plane) of a turbulent impinging 
jet. The jet temperature is 300 K. 
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Figure 10. Velocity distribution in the ground plane (x-z plane) of a turbulent impinging 
jet. The jet temperature is 1000 K. 
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